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ABSTRACT 

A  radar  target,  acting  as  a  scatterer  of  an  incident 
electromagnetic  wave,  can  be  considered  as  a  linear  time- 
invariant  system.  Previous  work  has  shown  ^hat  the  target's 
pole  locations  are  independent  of  the  incident  electromagnetic 
excitation,  including  incident  wave  shape,  aspect  and 
polarization.  This  thesis  develops  the  Kumaresan-Tuf ts  and 
Cadzow-Solomon  signal  processing  algorithms  into  computer 
routines  and  evaluates  their  pole  extraction  performance. 
Data  used  to  evaluate  the  extraction  algorithms  includes 
synthetic  and  integral  equation  generated  signals  with 
additive  noise,  in  addition  to  measurements  of  scattering  by 
scale  models  made  in  an  anechoic  chamber. 

J 


TABLE  OF  CONTENTS 


I .  INTRODUCTION  .  1 

A.  THE  PROBLEM . 2 

B .  BACKGROUND . 3 

C.  HISTORY  . 8 

II.  POLE  EXTRACTION  ALGORITHMS  .  9 

A.  PREVIOUS  WORK . 9 

1.  Direct  Minimization  .  9 

2.  Prony '  s  Method . 11 

B.  KUMARESAN-TUFTS  ALGORITHM  .  12 

1 .  Equations . 13 

2.  Singular  Value  Decomposition  ....  14 

3.  Bias  Compensation . 16 

4.  Kumaresan  and  Tufts  Compensation  .  .  16 

5.  Compensation  Based  on  Eigenvalue 

Shifting  Theorem  .  17 

6.  Performance . 19 

a.  Synthetically  Generated  Data  .  .  19 

1.  Noise  Performance  .  19 

b.  Thin  Wire  Integral  Equation 

Generated . 29 

c.  Scale  Model  Measurements  ....  38 

1 .  Wire  Targets . 50 

2.  Aircraft  Models . 50 

C.  CADZOW-SOLOMON  ALGORITHM  .  63 

1.  Applicability . 68 

2.  Equations . 68 

3.  Excess  Poles  and  Noise  Removal  ...  69 

4.  Singular  Value  Decomposition  ....  70 

5.  Bias  Compensation  in  the  Cadzow-Solomon 

Formulation . 70 

6.  Performance . 72 

a.  Synthetically  Generated  Data  .  .  72 

1.  Noise  Performance  .  72 

b.  Thin  Wire  Integral  Equation 

Generated  Data . 80 

c.  Scale  Models . 80 

1.  Wire  Targets . 92 

2.  Model  Aircraft . 92 


IV 


III.  SUMMARIES  AND  CONCLUSIONS 


108 


A.  KUMARESAN-TUFTS . 108 

f  B.  CADZOW- SOLOMON . Ill 

C.  CONCLUSIONS . 112 

APPENDIX  A.  THE  FUMARESAN-TUFTS  POLE  EXTRACTION 

ALGORITHM . 113 

APPENDIX  B.  THE  CADZOW-SOLOMON  POLE  EXTRACTION 

ALGORITHM . 125 

APPENDIX  C.  MATRIX  MULTIPLICATION  .  140 

APPENDIX  D.  GRAPHICS  ROUTINE  .  141 

LIST  OF  REFERENCES . 145 

INITIAL  DISTRIBUTION  LIST  .  147 


V 


LIST  OF  FIGURES 


Figure 

1. 

Signal  Containing  two  S-Plane  Poles, 

90.0  dB  SNR  . 

21 

Figure 

2  . 

Kumaresan-Tufts  Poles,  Synthetic  Data, 
90.0  dB  SNR  . 

22 

Figure 

3  . 

Kumaresan-Tufts  Poles,  Synthetic  Data, 
30.0  dB  SNR  . 

23 

Figure 

4  . 

Kumaresan-Tufts  Poles,  Synthetic  Data, 
20.0  dB  SNR  . 

24 

Figure 

5  . 

Kumaresan-Tufts  Poles,  Synthetic  Data, 
10.0  dB  SNR  . 

25 

Figure 

6  . 

Kumaresan-Tufts  Poles,  Synthetic  Data, 

7 . 0  dB  SNR  . 

26 

Figure 

7  . 

Signal  Containing  Two  S-Plane  Poles, 

7 . 0  dB  SNR  . 

2  7 

Figure 

8  . 

Kumaresan-Tufts  Pole  Extraction, 

7.0  dB  SNR  . 

28 

Figure 

9  . 

Double  Gaussian  Pulse  . 

30 

Figure 

10. 

Integral  Equation  Thin  Wire  Scattering, 
30  Degree  Aspect  . 

32 

Figure 

11 . 

Integral  Equation  Thin  Wire  Scattering, 
45  Degree  Aspect  . 

33 

F igure 

12  . 

Integral  Equation  Thin  Wire  Scattering, 
60  Degree  Aspect  . 

34 

Figure 

13  , 

Integral  Equation  Thin  Wire  Scattering, 
90  Degree  Aspect  . 

35 

Figure 

14  . 

Kumaresan-Tufts  Poles,  Noiseless  Thin 
Wire  Data  . 

36 

Figure 

15. 

Integral  Equation  Thin  Wire  Scattering, 
20.0  dB  SNR,  45  Degree  Aspect  .... 

37 

Figure 

16. 

Kumaresan-Tufts  Poles,  20.0  dB  SNR  .  . 

39 

VI 


Figure  17.  Integral  EguaLion  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

3  0  Degree  Aspect . 4  0 

Figure  18.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

4  5  Degree  Aspect . 41 

Figure  19.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

60  Degree  Aspect . 42 

Figure  20.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

90  Degree  Aspect . 4  3 

Figure  21.  Integral  Equation  Thin  Wire  Scattering, 

7.0  dB  SNR,  4  5  Degree  Aspect . 4  4 

Figure  22.  Kumaresan-Tufts  Poles,  7.0  dB  SNR  ...  45 


Figure  23.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

3  0  Degree  Aspect . 4  6 

Figure  24.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

4  5  Degree  Aspect . 4  7 

Figure  25.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

60  Degree  Aspect . 48 

Figure  26.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

90  Degree  Aspect . 49 

Figure  27.  Measured  Thin  Wire  Scattering, 

30  Degree  Aspect . 51 

Figure  28.  Measured  Thin  Wire  Scattering, 

45  Degree  Aspect . 52 

Figure  29.  Measured  Thin  Wire  Scattering, 

60  Degree  Aspect . 53 

Figure  30.  Measured  Thin  Wire  Scattering, 

90  Degree  Aspect . 54 

Figure  31.  Kumaresan-Tufts  Poles, 

Measured  Thin  Wire . 55 


VI 1 


f 


Figure  32.  Thin  Wire  Comparison,  Measured  vs. 

Integral  Equation  .  56 

Figure  33.  Target  1  Scattering,  30  Degrees 

from  Nose  on . 57 

Figure  34.  Target  1  Scattering,  Nose  on . 58 

Figure  35.  Target  2  Scattering,  30  Degrees 

from  Nose  on . 59 

Figure  36.  Target  2  Scattering,  Nose  on . 60 

Figure  37.  Kumaresan-Tuf ts  Poles,  Target  1  ...  .  61 

Figure  38.  Kumaresan-Tuf ts  Poles,  Target  1  .  .  .  .  63 

Figure  39.  Kumaresan-Tuf ts  Poles,  Target  1  .  .  .  .  64 

Figure  40.  Kumaresan-Tuf ts  Poles,  Target  2  .  .  .  .  65 

Figure  41.  Kumaresan-Tuf ts  Poles,  Target  2  ...  .  66 

Figure  42.  Kum.arecan-Tuf ts  Poles,  Target  2  .  .  .  .  67 

Figure  43.  Signal  Containing  Two  S-Plane  Poles, 

90.0  dB  SNR . 73 

Figure  44.  Cadzow-Solomon  Poles,  Synthetic  Data, 

0  dP  SNP .  ...  74 

Figure  45.  Cadzow-Solomon  Poles,  Synthetic  Data, 

30.0  dB  SNR . 75 

Figure  46.  Cadzow-Solomon  Poles,  Synthetic  Data, 

20.0  dB  SNR .  . 

Figure  47.  Cadzow-Solomon  Poles,  Synthetic  Data, 

10.0  dB  SNR . 77 

Figure  48.  Cadzow-Solomon  Poles,  Synthetic  Data, 

7.0  dB  SNR . 78 

Figure  49.  Cadzow-Solomon  Poles,  Noiseless 

Thin  Wire  Data . 81 

Figure  50.  Cadzow-Solomon  Poles,  20.0  dB  SNR  ...  82 

Figure  51.  Integral  Equation  Thin  Wire  Comparison, 

Noiseless  vs.  20.0  dB  SNR, 

30  Degree  Aspect . 83 

viii 


Figure  52.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

45  Degree  Aspect . 84 

Figure  53.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

60  Degree  Aspect . 85 

Figure  54.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  20.0  dB  SNR, 

90  Degree  Aspect . 86 

Figure  55.  Cadzow-Solomon  Poles,  7.0  d5  SNR  ...  87 

Figure  56.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

30  Degree  Aspect . 88 

Figure  57.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

45  Degree  Aspect . 89 

Figure  58.  Integral  Equation  Thin  Wire  Comparison. 
Noiseless  vs.  7.0  dB  SNR, 

oO  Degree  Aspect . 90 

Figure  59.  Integral  Equation  Thin  Wire  Comparison, 
Noiseless  vs.  7.0  dB  SNR, 

90  Degree  Aspect . 91 

Figure  60.  Cadzow-Solomon  Poles, 

Measured  Thin  Wire . 97 

Figure  61.  Thin  Wire  Comparison, 

Measured  vs.  Integral  Eauation  ....  94 

Figure  62.  Cadzow-Solomon  Poles  Target  1, 

Three  Aspects . 95 

Figure  63.  Cadzow-Solomon  Poles  Target  1, 

Three  Aspects . 96 

Figure  64.  Cadzow-Solomon  Poles  Target  1, 

All  Six  Aspects . 97 

Figure  65.  Cadzow-Solomon  Poles  Target  2, 

Three  Aspects . 98 

Figure  66.  Cadzow-Solomon  Poles  Target  2, 

Three  Aspects . 99 


IX 


Fi yure  67  . 


Cadzow-Solomon  Poles  Target  2, 
All  Six  Aspects  . 


100 


I 

r 

I 


Figure  68.  Pole  Comparisons,  Target  1, 

All  Six  Aspects . 101 

Figure  69.  Pole  Comparisons,  Target  2, 

All  Six  Aspects . 102 

Figure  70.  Target  3  Scattering,  Nose-on  .  104 

Figure  71.  Target  4  Scattering,  Nose-on  .  105 

Figure  72.  Cadzow-Solomon  Pole  Comparisons, 

4  Targets,  Nose  on . 106 


I 

i 

1 

r 


I.  INTRODUCTION 


A  radar  target,  acting  as  a  scatterer  of  a  specified 
incident  electromagnetic  wave,  can  be  considered  as  a  single 
input,  single  output,  linear  time-invariant  (LTD  system  for 
a  fixed  field  observation  point.  The  target  can  thus  be 
considered  as  a  transfer  function  with  poles  and  zeros.  Baum 
demonstrated  at  the  Air  Force  Weapons  Laboratory  that  a 
target's  induced  current  response  ro  an  incident  electro¬ 
magnetic  wave  has  identifiable  poles  determined  by  the 
composition  and  structural  geometry  of  the  target  [1].  In 
1974,  Moffatt  and  Mains  proposed  that  the  target's  scattered 
field  pole  locations  are  independent  of  the  incident 
e]  ectrom.agnetic  excitation,  including  aspect  and  polarization 
[2].  Morgan  has  proven  theoretically  that,  for  the  case  of 
a  conducting  target,  the  scattering  response  contains  complex 
natural  resonances  which  are  independent  of  the  incident 
electromagnetic  excitation  [3].  By  determining  the  poles  of 
a  target's  response,  aspect  independent  target  identification 
can  be  accomplished  through  the  use  of  electromagnetic  natural 
resonances . 

Although  the  concept  of  radar  target  identification 
through  the  use  of  natural  resonances  was  first  proposed  in 
1974  by  Mains  and  Moffatt  [2],  only  recently  have  signal 
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processing  techniques  been  applied  to  locate  the  poles  in  a 
radar  target's  response  in  the  presence  of  noise.  Kumaresan- 
Tufts  [4]  and  Cadzow-Solomon  [5]  have  each  developed 
algorithms  which  have  proven  successful  in  the  presence  of 
noise.  This  thesis  develops  computer  routines  based  upon 
these  two  algorithms  and  examines  their  respective  performance 
and  appropriateness  using  a  variety  of  scattering  data. 

A.  THE  PROBLEM 

Since  the  performance  of  signal  processing  methods  varies 
under  different  conditions,  a  system  employed  to  identify 
targets  would  possibly  reach  a  decision  based  on  the  combined 
output  of  several  signal  processing  methods.  For  example,  the 
Kumaresan-Tuf ts  and  Cadzow-Solomon  methods  could  be  used  to 
extract  poles  from  the  response  of  scale-  model  targets.  The 
information  so  ga+-hered  could  be  used  to  build  a  data  base  for 
comparison  with  data  similarly  obtained  in  actual  field  use. 
The  results  of  this  system  would  serve  as  one  input  to  a 
larger  system.  Other  methods  would  provide  input  to  the 
system,  such  as  the  K-pulse  method  of  Kennaugh  [6]  and  the 
annihilation  filter  used  by  Dunavin  [7],  Morgan  and  Dunavin 
[8]  and  Chen  [9].  As  the  name  suggests,  an  annihilation 
filter  annihilates  the  target's  poles.  A  system  using  the 
annihilation  filter  concept  would  contain  many  such  filters, 
each  previously  designed  to  cancel  the  poles  of  a  specific 
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known  target.  In  actual  field  use,  a  radar  target's  response 
would  be  input  into  each  of  the  filters,  and  the  target 
selected  would  be  that  matching  the  filter  whose  output 
exhibits  the  lowest  signal  energy. 

A  system  used  to  identify  radar  targets  would  require  the 
following  concept  of  employment.  First,  information  required 
by  each  of  the  sub-systems  would  be  obtained  for  every  target 
class  of  concern.  In  actual  field  use,  this  information  would 
be  compared  against  actual  radar  target  responses.  The  system 
would  then  determine  the  identity  of  the  target  based  on  the 
input  from  each  of  its  sub-systems. 

B .  BACKGROUND 

Consider  a  perfectly  conducting  target  illuminated  by  an 
electromagnetic  field.  The  current  induced  on  the  surface  of 
this  target  at  a  given  point  must  satisfy  the  magnetic  field 
integral  equation  (MFIE) , [10] 

j(r,i)=2nxH‘(r,0  +  f  J  K  (?,  r  ,  f)  J  ( r I  )  dS  dj 

Spv 

where  h  is  an  outward  unit  vector  normal  to  the  surface  of 
the  object,  J  is  the  surface  current  density,  h',  is  the 
incident  magnetic  field,  and  K  is  a  Green's  function  dyadic. 
The  entire  equation  is  most  easily  understood  as  the  sum  of 
driven  currents  and  "feedback"  currents  corresponding  to  the 
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cross-product  term  and  surface  integral  term  respectively. 
The  term  driven  by  the  magnetic  f ield , 2nxH , forms  the  physical 
optics  portion  of  the  total  current.  Physical  optics 
describes  the  cross-product  term  as  the  induced  current 
without  interaction  with  the  rest  of  the  body.  The  Green's 
function  kernel  describes  the  current  at  a  point  on  the  object 
due  to  the  feedback  of  currents  from  every  other  point  on  the 
object,  as  previously  illuminated  by  the  incident  field.  The 
current  at  each  point  is  then  summed  over  the  surface  of  the 
object.  Note  that  the  surface  integral  term  is  of  principal- 
value  type;  the  integral  excludes  the  point  f=f  . 

Once  the  incident  magnetic  field  is  no  longer  present, 
the  solutions  of  (1)  are  considered  the  natural  modes  of  the 
object.  These  natural  modes  are  of  the  form,  J^exp(s^)  .  The 
natural  resonance  frequencies  are  of  the  form, 

(2) 

where  is  the  damping  rate  in  Nepers/sec  and  0)^  is  the 
frequency  in  radians/sec.  The  natural  resonances  of  (2)  are 
functions  of  the  structural  geometry  of  the  object  and  are 
independent  of  the  incident  magnetic  field.  To  understand 
how  these  natural  resonances  are  unique  to  the  geometry  and 
composition  of  the  object,  consider  a  set  of  points  on  the 
object  previously  illuminated  by  the  incident  field,  so  that 
h'=0.  The  current  at  a  given  point  in  the  set  is  due  to  the 


infinite  number  of  feedback  currents  from  every  other  point 


in  the  set.  Recall  that  these  feedbacks  are  described  by  the 
Green's  function  kernel  in  the  integral  term  of  (1).  Since 
the  set  of  points  previously  illuminated  is  physically  located 
on  the  same  object,  the  infinite  number  of  paths  that  connect 
a  point  with  all  other  points  in  the  set  is  the  same  for  all 
points  in  the  set.  The  infinite  number  of  paths  are  unique 
to  the  structural  geometry  of  the  object  and  correspond 
exactly  to  the  infinite  number  of  paths  taken  by  currents 
which  feedback  to  a  given  point  via  the  Green's  function 
kernel.  Finally,  the  composition  of  the  target  determines  the 
surface  current  density  on  the  object.  Although  an  infinite 
number  of  resonances  exists  in  any  object,  only  a  limited 
number  of  these  will  be  measurably  excited  by  an  incident 
field  of  finite  bandwidth.  These  resonances  described  in  (2) 
appear  as  complex  conjugate  pairs  in  the  left-half  portion  of 
the  s-plane. 

In  the  far-field,  the  back-scattered  response  of  a  target 
to  an  incident  plane  wave  is  of  the  form 

t)  (r  ,  t-|f-r  |/c)  dS'  (3 ) 

where  c  is  the  speed  of  light  and  p  is  the  unit  vector  whose 
direction  matches  that  of  the  plane  wave's  propagation. 
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Equation  (3)  is  the  result  of  integrating  the  current  at 
each  point  on  the  target  surface  for  a  fixed  point  in  the  far- 
field.  Recall  that  the  current  at  each  point  on  the  target 
is  defined  by  (1).  Thus,  the  back-scattered  far-field  can 
be  obtained  by  substituting  (1)  into  (3): 

H  (-rp,t)=u(t-r/c)  f  H  (-rp,t)+  X  H  J-rp,  t)  exp  (s  t) 

I  n=-“ 

^  n?:o 

The  currents  in  (1)  produce  the  field  in  (4).  In  fact,  each 
term  in  (4)  corresponds  to  the  term  in  (1)  which  produced  it. 
Specifically,  the  first  term  in  (4)  describes  the  physical 
optics  scattered  field  generated  by  the  2nxH'  current  which, 
of  course,  is  the  first  term  in  (1).  Similarly,  the  second 
term  in  (4)  is  produced  by  the  source-free  currents  defined 
by  the  second  term  in  (1) .  Like  the  current  described  in  (1) , 
the  field  in  (4)  is  the  sum  of  two  terms,  a  driven  term  and 
a  term  containing  feedbacks. 

The  results  of  (4)  can  also  be  seen  as  two  forms  of  the 
Singularity  Expansion  Method  (SEM)  developed  by  Baum  [1]  .  As 
shown  by  Morgan  [10] ,  during  the  early-time  portion  of  the 
target's  response,  the  scattered  field  is  composed  of  the 
physical  optics  scattered  field  and  a  "Class  2"  form  of  the 
SEM  expansion.  The  class  2  SEM  expansion  corresponds  to  the 
second  term  of  (4),  wherein  the  coefficients  are  time- 
varying  as  the  wave  passes  over  the  target,  since  the  currents 
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producing  this  portion  of  the  field  are  integrated  over  a 
time-varying  surface  area.  At  the  instant  the  wave  passes  the 
last  point  of  the  target,  the  physical  optics  field  vanishes 
and  the  remaining  term  in  (4)  is  produced  by  constant 
coefficients  .  The  coefficients  are  constant  at  this 
instant  since  the  surface  area  in  the  integral  in  (3)  is  now 
constant.  This  instant  also  marks  the  transition  of  (4)  from 
a  "class  2"  SEM  expansion  of  time-varying  coefficients  to  a 
"class  1"  SEM  expansion  of  constant  coefficients.  The 
scattered  field  due  to  a  plane  wave  is  therefore  composed  of 
a  physical  optics  term  and  a  class  2  SEM  expansion  in  the 
early-time,  and  a  simple  class  1  expansion  in  the  late-time. 

Actual  measurement  of  the  scattered  far-zone  field  would 
be  greatly  aided  by  knowledge  of  the  transition  time  of  the 
field  from  early  time  to  late  time.  From  [10],  this 
transition  for  a  monostatic  radar  would  occur  at  At=T+2 (D+d) /c 
seconds  after  radar  turn-on.  Here,  T  is  the  pulse  duration, 
D  is  the  target's  dimension  along  the  direction  of  wave 
propagation,  d  is  the  distance  between  the  target  and  the 
measurement  point  and  c  is  the  speed  of  light. 

The  discussion  presented  in  this  section  was  extracted 
from  work  done  by  Morgan  in  [10] .  The  reader  is  referred  to 
this  work  for  a  more  detailed  treatment  of  the  material  in 
this  section. 


I 


7 


HISTORY 


C  . 

The  results  of  the  previous  section  form  the  basis  for  the 
hypothesis  that  the  natural  resonances  found  in  the  scattering 
response  of  a  target  to  an  incident  electromagnetic  wave  are 
unique  to  that  target.  Additionally,  only  a  finite  set  of 
these  natural  resonances  are  measurably  excited  by  a  wave  of 
finite  bandwidth.  In  1974,  Moffatt  and  Mains  proposed  that 
the  extraction  of  resonances  from  a  target's  response  to 
electromagnetic  excitation  could  be  used  for  target 
identification.  This  work  related  to  earlier  work  in  1965, 
when  Kennaugh  and  Moffatt  first  developed  the  concept  of  a 
radar  target  as  a  linear  time  invariant  system.  Poles  in  the 
z-plane  are  directly  related  to  the  natural  resonances  of  a 
target 

z„-e  " 

"  (5) 

where  is  given  by  (2)  and  At  is  the  sampling  interval  in 
seconds.  Hence,  pole  extraction  involves  resonance 

identification.  The  use  of  pole  extraction  algorithms  is 
discussed  in  the  next  chapter. 
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II.  POLE  EXTRACTION  ALGORITHMS 


The  use  of  pole  extraction  algorithms  to  identify  radar 
targets  is  discussed  in  this  chapter.  A  brief  discussion  of 
two  methods  precedes  the  in-depth  evaluation  of  the  Kumaresan- 
Tufts  and  Cadzow-Solomon  algorithms.  The  evaluation  of  the 
latter  two  algorithms  occurs  in  two  stages.  First,  each 
algorithm  will  be  evaluated  in  its  ability  to  extract  poles 
from  data  with  known  poles.  Some  of  the  data  processed  was 
generated  at  various  signal  to  noise  ratios  by  a  computer 
program  written  by  Morgan  [11] .  Additional  data  was  produced 
by  Morgan's  time-domain  thin  wire  integral  equation  computer 
program  [12].  In  the  second  stage,  a  side  by  side  comparison 
is  made  of  poles  extracted  by  each  method  using  transient 
scattering  measurements  for  a  thin  wire  and  for  various  model 
aircraft.  Comparisons  between  the  two  methods  are  made  as  the 
aspect  of  the  aircraft  is  varied. 

A.  PREVIOUS  WORK 

1.  Direct  Minimization 

The  most  direct  way  to  determine  the  natural 
resonances  in  a  target’s  response  is  to  minimize  the  mean- 
square  error  between  the  modeled  signal  and  the  received 
signal.  In  [10],  Morgan  determined  that  the  late-time  target 
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response  to  a  radar  could  be  represented  as  a  sum  of  damped 
sinusoids  given  by 


“  0  I 

y(t)  =  2  A,e  ‘  cos(w,t+ei)  (6) 

1  =  1 

The  frequency,  u,  ,  and  damping  rate,  „  ,  are  the  same 
parameters  found  in  the  natural  resonance  defined  in  (2). 
Phase,  6|  ,  and  amplitude.  A,,  are  the  remaining  parameters. 
The  representation  in  (6)  is  the  sum  of  an  infinite  number  of 
resonances.  The  sampled  response  to  an  incident  wave  of 
finite  bandwidth  can  be  modeled  as 

N  0 ,  na  i 

y(nat)  =  y=SA.e  cos  (w,nat+0,)  (7) 

”  1  =  1 

where  At  is  the  sampling  interval  in  seconds.  The  four 
parameters  of  (7)  must  be  adjusted  to  minimize  the  sampled 
mean-square  error  signal 

®n=(y,-yn)'  (8) 

between  the  actual  discrete  sampled  received  signal  y  and  the 

n 

modeled  signal  y^ .  The  processing  required  in  this 
minimization  problem  is  both  inefficient  and  highly  non¬ 
linear.  Nevertheless,  Chong  used  this  method  to  process 
mathematically-generated  data  down  to  15.0  dB  signal-to-noise 
( SNR)  ratio  [13]  . 
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2.  Prony's  Method 

As  in  direct  minimization,  Prony's  approach  to 
resonance  classification  focuses  on  the  late-time  portion  of 
a  radar  target's  response.  However,  linear  processing  and 
root  solving  are  used.  The  late-time  response  is  modeled  as 
the  output  of  an  LTI  system  of  order  k^.-  Each  signal  received 
at  some  discrete  sample,  n,  is  considered  to  be  the  weighted 
sum  of  Kp  previous  signals.  Thus,  the  finite  term 
approximation  of  the  received  late-time  signal,  y  ,  is  defined 


The  z-transform  of  (9)  is 

The  roots  of  this  polynomial  in  z  are  the  poles  of  the  system 
model.  Therefore,  the  key  to  extracting  the  poles  in  the 
system's  response  lies  in  solving  for  the  coefficients  b  of 


A  set  of  Kp+M  received  signals  in  M  equations  (9)  can 
be  arranged  in  matrix  form  as 


■  n  •••  'V'  fb,  -|  r  y 

•  •  •  ,  ^ 
0  •  •  =s  • 
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In  Prony's  original  method,  the  data  matrix  is 
exactly  determined,  and  the  coefficient  vector,  b,  is  solved 
using  linear  computations.  In  the  presence  of  noise,  Prony 
overdetermines  the  data  matrix  by  setting  M>Kp  and  solves  for 
the  coefficient  vector  by  obtaining  the  least-squares  solution 
to  the  system  of  equations. 

The  Prony  method  has  two  major  problems.  First,  the 
poles  obtained  by  the  least  squares  solution  to  the 
overdetermined  matrix  may  be  strongly  perturbed  by  noise  [14] , 
since  noise  does  not  satisfy  the  causal  model  of  the  system. 
Second,  the  order  of  the  system  is  generally  not  known  a 
priori.  When  the  estimated  order  is  greater  than  the  actual 
order,  poles  due  to  noise  are  generated.  Prony's  method 
offers  no  technique  for  distinguishing  between  the  signal 
poles  and  the  extra  poles  caused  by  overestimation  of  the 
system's  order.  If  the  estimated  system  order  is  less  than 
the  actual  order,  actual  poles  are  lost  and  the  remaining 
poles  are  perturbed  from  their  true  positions. 

B.  KUMARESAN-TUFTS  ALGORITHM 

The  Kumaresan  and  Tufts  pole  extraction  algorithm  was 
developed  by  adapting  Prony's  method  to  reduce  the  problems 
addressed  in  the  preceding  section.  The  Kumaresan-Tuf ts 
algorithm  modifies  the  least-squares  Prony  method  in  three 
ways  : 
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1.  Processed  signals  are  arranged  in  a  data  matrix 
based  on  a  non-casual  model  of  the  system. 

2.  The  model  of  the  system  is  deliberately 
overestimated . 

3.  The  system  of  equations  determined  by  the  above 
two  criteria  is  solved  by  using  singular  value 
decomposition  (SVD). 

Kumaresan  demonstrates  in  [15]  that  the  use  of  singular 
value  decomposition  tends  to  force  the  extra  poles  of  the 
excess-order  system  inside  the  unit  circle,  while  the  non- 
causal  arrangement  of  the  signals  tends  to  force  the  signal 
poles  outside  the  unit  circle.  The  excess  order  of  the  system 
model  reduces  the  effects  of  noise  on  the  actual  poles.  Since 
the  noise  is  stationary  and  stable,  it  looks  the  same  in 
forward  and  backward  time. 

1 .  Equations 

Recall  that  in  (9),  Prony ’ s  technique  defines  the 
received  late-time  signal  as  the  weighted  sum  of  previous 
signals,  where  is  presumed  to  be  the  order  of  the  system. 
Kumaresan  models  the  same  late-time  signal  as  the  weighted  sum 
of  Kp  future  signals,  where  Kp  is  greater  than  the  estimated 
order  of  the  system.  This  non-casual  model  is  given  by 


(12) 
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As  in  Prony's  method,  the  coefficients  b',  are  coefficients  of 
a  polynomial  in  z  that  models  the  system's  late-time  response. 
Two  simple  manipulations  of  either  data  matrix  leads  to  the 
relationship  between  the  coefficients  of  the  Prony  model  and 
the  prediction  coefficients  of  the  Kumaresan-Tuf ts  model. 
With  b  =-i ,  a  prediction  coefficient  is  related  to  an 
autoregressive  coefficient  by 


From  the  above  relationship,  it  can  be  shown  that  the  complex 
pole  pairs  of  the  causal  model  are  merely  conjugate 
reflections  across  the  unit  circle  of  the  pole  pairs  in  the 
non-causal  model. 

2.  Singular  Value  Decomposition 

The  non-causal  arrangement  of  late-time  signals  in  a 
set  of  system  equations,  and  subsequent  processing  through 
singular  value  decomposition,  ccm.bine  to  separate  the  signal 
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and  noise  into  orthogonal  spaces.  As  discussed  in  the 
preceding  paragraph,  poles  of  the  non-causal  model  are 
reflected  outside  the  unit  circle.  Kumaresan  demonstrates  in 
[15]  that  the  extra  poles  of  the  excess-order  system  can  be 
forced  inside  the  unit  circle  through  the  use  of  SVD . 

Singular  value  decomposition  factors  the  -  MXKj,'  data 
matrix  D.^.  into  the  product  of  the  matrices: 


The  columns  of  U  (MXM)  are  eigenvectors  of 

columns  of  V  (K^XKp)  are  eigenvectors  of  .  If  r  is  the 

rank  of  the  data  matrix,  ,  the  diagonal  matrix  s  (MXK.,) 
contains  r  singular  values  which  are  the  square  roots  of  the 
nonzero  eigenvalues  of  both  and  .  By  rearranging 

the  three  matrices  in  the  product,  the  pseudoinverse  of  ^y.  can 
be  obtained  as 


Dy=Vi*U'^ 


(17) 


where  x*  is  a  (K^XM)  matrix  whose  singular  values  on  the 
diagonal  are  the  reciprocals  of  those  in  the  I  matrix. 
Finally,  the  coefficient  vector  b*,  of  minimum  Euclidian  norm, 
is  given  by 
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The  coefficient  vector  b'"  so  obtained  is  the  minimum  length 
least-squares  solution  to  (14).  In  other  words,  b*  is  the 
best  possible  solution  to  (14)  .  In  the  case  of  noiseless 
data,  the  extraneous  poles  generated  by  the  excess-order  model 
will  always  be  inside  the  unit  circle  when  b*  is  used.  This 
result  is  generally  true  for  noisy  data. 

3.  Bias  Compensation 

Kumaresan  and  Tufts  [4]  observed  that  the  addition  of 
noise  perturbed  the  singular  values  of  the  z  matrix  of  (16) . 
If  the  perturbation  of  these  singular  values  is  not 
compensated,  both  the  signal  poles  and  extraneous  poles  are 
biased  towards  the  unit  circle.  Kumaresan  and  Tufts  used  a 
compensation  method  which  reduced  the  bias  in  their  wor)c,  but 
did  not  derive  an  analytical  justification.  In  [16],  Norton 
derived  a  more  valid  bias  compensation  method  based  on  the 
eigenvalue  shifting  theorem. 

4.  Kumaresan  and  Tufts  Compensation 

If  the  actual  order  of  the  system  is  then  the 
first  Kj,  singular  values  of  the  I  matrix  in  (16)  are  non¬ 
zero.  The  remaining  Kp-Kj,  singular  values  are  considered 
noise  singular  values  and  are  zero  in  the  case  of  noiseless 
data.  The  addition  of  noise  perturbs  the  first  Kp  signal 
singular  values  and  increases  the  noise  to  some  non-zero 
value.  Kumaresan  and  Tufts  compensated  for  this  increase  in 
the  singular  values  due  to  the  noise  by  subtracting  the 
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average  of  the  noise  singular  values  from  the  signal  singular 
values.  The  noise  singular  values  were  then  set  to  zero. 

5.  Compensation  Based  on  Eigenvalue  Shifting  Theorem 

As  described  in  the  previous  section,  the  singular 
values  of  the  matrix  are  the  square  roots  of  the 

eigenvalues  of  DyDy.  and  D^Dy  .•  Assume  the  noisy  data  matrix 
can  be  represented  by  d^=S+N'/  where  N  is  composed  of  the  wide- 
sense  stationary  white  noise  process  v,,  given  by 


The  expected  value  of  D^Dy  can  be  obtained  by 

DyD^=E[  (S  +  N)  (S  +  N)^]=ECss’']+E[SNn+ELNS’']+E[NN’']  (20) 


Since  S  is  deterministic,  ECss’^  ]=SS'^ .  Assuming  the  noise  is 
zero  mean,  the  two  cross  products  are  zero.  Because  we  assume 
the  noise  is  wide-sense  stationary  and  white,  E[NN^]=a^I, 
where  is  the  noise  variance  and  I  is  the  identity  matrix. 
The  expected  value  of  !  D  thus  becomes 

y  y- 


E[DyD^]=SS^+o2l 


(21) 
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Similarly,  the  expected  value  of  d^D  «  the  other  source  of 

y  y 

singular  values,  is 

^22) 

The  assumption  in  the  results  of  (21)  and  (22)  is  that  the 
diagonals  of  E  [N^N] =E [NN^ ]  equals  the  noise  variance  cl  .• 
Equations  (21)  and  (22)  show  that  in  the  mean,  the  squares  of 
the  singular  values  of  are  increased  by  the  noise  variance. 

The  results  lead  to  the  method  of  eigenvalue 

compensation  recommended  by  Norton  in  [16].  Recall  from  (16) 
that  the  eigenvalues  of  Dy  are  on  the  diagonal  of  the  Z 
matrix  returned  by  the  singular  value  decomposition  of  Dy  • 
If  Kp'  is  the  actual  order  of  the  system,  and  Kp  is  the 

estimated  order  of  the  system  then  the  remaining  Kp-Kp' 
singular  values  of  the  Z  matrix  can  be  squared  and  averaged 
to  obtain  an  estimate  of  the  noise  variance,  •  These  noise 

singular  values  can  then  be  set  to  zero.  The  first  Kp 

singular  values  of  the  Z  matrix  are  then  squared  and  reduced 
by  subtracting  the  estimate  of  the  noise  variance.  The  square 
root  of  the  difference  becomes  the  new  first  Kp  singular 
values  of  the  compensated  Z  matrix.  Calculations  according 
to  (17)  and  (18)  can  then  be  carried  out  in  a  normal  manner 

to  obtain  poles  in  the  presence  of  the  noise.  Eigenvalue 

compensation  requires  an  estimate  of  the  actual  order  of  the 
system.  Methods  to  obtain  this  estimate  are  discussed  in 
Chapter  III. 
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6.  Performance 

The  Kumaresan-Tuf ts  algorithm  was  programmed  in 
Fortran  and  tested  on  various  types  of  data.  The  program 
appears  in  Appendix  A. 

a.  Synthetically  Generated  Data 

The  starting  point  for  evaluating  the  performance 
of  the  Kumaresan-Tuf ts  algorithm  was  with  synthetically 
generated  data  of  the  form  given  by  (8)  and  shown  here  again 
for  convenience 


,  O.nat 

‘  cos ((i>|mt+9, )  (8) 

Again,  A, ,  a, ,  c, , 0^ ,  are  the  amplitude,  damping  rate,  frequency 
and  phase  of  a  set  of  N  damped  sinusoids.  Noisy  data  was 
created  by  adding  stationary  white  noise. 

1.  Noise  Performance 

The  algorithm  was  evaluated  at  various  SNR's, 
ranging  from  90.0  dB  to  7.0  dB.  These  SNR's  are  ratios  of 
signal  energy  to  noise  energy  rather  than  the  ratio  of  signal- 
to-noise  power.  Synthetic  data  so  generated  more  closely 
resembles  the  exponential  decay  of  signal  power  typical  in 
actual  radar  measurements. 
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Figure  1  shows  the  signal  produced  by  two  s- 
plane  poles  at  90.0  dB.  Figures  2  through  6  depict  the  poles 
extracted  from  this  signal  at  SNR's  ranging  from  90.0  dB  to 
7.0  dB.  Obtained  poles  are  shown  at  their  positions  within 
the  upper  right  hand  quadrant  of  the  unit  circle  in  the  z- 
plane.  Not  shown  are  conjugates  of  each  pole  which  are 
located  below  the  real  axis  outside  the  figure  boundaries. 

Figures  2  through  6  demonstrate  outstanding 
performance  on  noisy  data,  even  at  SNR's  of  7.0  dB .  The 
scaling  needed  to  show  a  discernible  difference  between 
results  obtained  at  30.0  dB  and  7.0  dB  would  necessarily 
exclude  one  of  the  poles  from  the  enlarged  figure.  The 
average  distance  of  the  trial  poles  obtained  in  the  7.0  dB 
SNR  signal  from  the  true  poles  is  on  the  order  of  10"^.  This 
magnitude  corresponds  to  that  of  the  average  estimate  of  the 
noise  variance  obtained  in  successive  trials  with  this  signal. 
The  correlation  between  the  distance  of  trial  poles  from  true 
poles  and  the  noise  variance  estimate  was  consistently 
observed  with  each  of  the  different  signal-to-noise  ratios 
used.  Figure  7  depicts  the  signal  of  Figure  1  severely 
corrupted  by  noise  having  7.0  dB  SNR. 

As  discussed  previously,  the  signal-to-noise 
ratio  used  in  the  synthetically  generated  data  is  the  ratio 
of  energy.  Figure  8  depicts  the  results  of  pole  extraction 
from  the  signal  shown  in  Figure  7,  but  with  a  late-time 
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Figure  2.  Kumaresan-Tufts  Poles,  Synthetic  Data,  90.0  dB  SNR 
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beginning  ten  nanoseconds  later.  Since  the  SNR  is  calculated 
over  twenty  nanoseconds  for  both  signals,  the  signal  power  at 
some  later  time  will  clearly  be  less  than  the  power  ten 
nanoseconds  earlier.  The  results  in  Figure  8  show  complete 
breakdown  of  the  algorithm's  ability  to  extract  poles.  The 
trial  poles  shown  are  the  poles  closest  to  the  true  poles,  and 
yet  they  are  located  at  positions  whose  reflections  are  inside 
the  unit  circle  where  noise  poles  are  typically  located. 

The  preceding  results  show  outstanding 
accuracy  for  full-length  noisy  data  but  a  complete  breakdown 
of  the  algorithm  for  the  same  signal  with  a  later  transition 
to  late-time.  These  initial  observations  are  supported  by 
similar  findings  presented  in  this  thesis. 

b.  Thin  Wire  Integral  Equation  Generated  Data 

For  simple  objects  such  as  a  thin  wire,  the  radar 
response  of  that  object  can  be  computed  by  establishing 
boundary  conditions  on  the  object  and  numerically  solving  the 
integral  equations  that  describe  the  surface  current.  Recall 
the  magnetic  field  integral  equation  given  by  (1). 
Simulations  produced  by  Morgan’s  time-domain  thin  wire 
integral  equation  computer  program  [12]  were  used  to  evaluate 
the  pole  extraction  algorithm.  The  excitation  waveform  used 
is  the  double  Gaussian  pulse  depicted  in  Figure  9.  This  pulse 
is  a  wide  Gaussian  pulse  with  a  ten  percent  width  of  0.3 
nanoseconds  subtracted  from  a  narrow  Gaussian  pulse  with  a  ten 
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Figure  9.  Double  Gaussian  Pulse 
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nanoseconds  subtracted  from  a  narrow  Gaussian  pulse  with  a  ten 
percent  width  of  0.15  nanoseconds. 

Figures  10  through  13  depict  back  scattering 
response  of  a  0.1  meter  length  thin-wire,  having  a  radius  of 
0.00118  meter,  computed  at  various  incident  aspects,  ranging 
from  thirty  degrees  to  ninety  degrees.  The  laboratory 
arrangement  for  actual  measurements  simulated  by  Morgan's 
program  is  described  in  [17]  .  Ninety  degrees  represents  a 
broadside  aspect,  while  thirty  degrees  represents  the  incident 
plane  wave  having  nearly  grazing  incidence  on  the  wire.  The 
poles  extracted  at  each  of  the  four  aspect  angles  are  plotted 
in  Figure  14.  In  this  figure,  and  those  that  follow  which 
depict  extracted  poles,  the  signal  poles  lie  in  or  on  the  unit 
circle,  and  the  noise  poles  lie  outside. 

The  results  obtained  with  this  rigorous  numerical 
computation  demonstrate  the  aspect  independence  of  the 
extracted  poles  using  the  Kumaresan-Tuf ts  method.  Note  that 
only  halt  of  the  poles  were  obtained  for  broadside 
illumination;  two  even-numbered  poles  can  easily  be  seen 
outside  the  unit  circle.  This  results  because  of  the  physical 
symmetry  of  both  the  wire  and  the  incident  field,  thus 
precluding  excitation  of  odd-symmetric  modal  currents  and 
their  associated  natural  resonances. 

Figure  15  exemplifies  the  computed  back-scattering 
response  of  the  0.1  meter  thin  wire  corrupted  artificially 
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Figure  13.  Integral  Equation  Thin  Wire  Scattering,  90  Degree 
Aspect 
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Figure  15.  Integral  Equation  Thin  Wire  Scattering,  20.0  dB 
SNR,  45  Degree  Aspect 


with  noise  at  a  20.0  dB  SNR.  Figure  16  shows  the  poles 
extracted  at  each  of  the  four  angles  of  incidence  used 
pieviously  in  Figure  14.  Poles  of  Figure  14  at  90°  are  now 
missing  in  Figure  16,  and  only  the  first  three  low  frequency 
poles  are  tightly  grouped.  The  loss  of  high  frequency  poles 
is  expected  because  these  have  the  highest  damping  and  thus 
lose  their  energy  at  the  fastest  rate.  Further  comparison 
between  results  computed  at  20.0  dB  SNR  and  infinite  SNR  are 
offered,  angle  by  angle,  in  Figures  17  through  20. 

One  additional  test  of  the  computed  thin  wire 
scattering  was  conducted  at  a  7.0  dB  SNR.  The  corrupted 
waveforms  are  exemplified  by  Figure  21;  the  extracted  poles 
are  shown  in  Figure  22.  The  number  of  poles  obtained  has 
decreased  with  respect  to  the  number  obtained  at  20.0  dB  SNR. 
The  grouping  of  the  clusters  has  also  expanded.  Angle  by 
angle  comparisons  are  again  offered  in  Figures  23  through  26. 
c.  Scale  Model  Measurements 

The  transient  scattering  measurements  of  scale 
models  used  for  evaluation  in  this  section  were  made  by  Walsh 
using  the  anechoic  chamber  of  the  Transient  Electromagnetic 
Scattering  Laboratory  at  the  Naval  Postgraduate  School.  The 
entire  measurement  process  and  laboratory  setup  are  described 
in  detail  in  [17] . 
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Figure  16.  Kumaresan-Tuf ts  Poles,  20. C  dB  SNR 
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Figure  17.  Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  30  Degree  Aspect 
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Integral  Equation  Thin  Wire  comparison,  Noiseless 
vs.  20.0  dB  SNR,  90  Degree  Aspect 
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Figure  21.  Integral  Equation  Thin  Wire  Scattering,  7.0  dB 
SNR,  4“^  Degree  Aspect 
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1.  wire  Targets 

The  thin  wire  measurements  were  obtained  from 
the  scattering  response  of  a  0.1  meter  length  thin  wire  having 
radius  0.00118  meter.  Recall  that  these  are  the  same 
dimensions  as  the  wire  whose  computed  response  was  processed 
in  the  previous  section.  The  measurements  at  each  of  four 
incident  aspects  are  shown  in  Figures  27  through  30. 

The  poles  extracted  from  the  four  measurements 
are  depicted  in  Figure  31.  As  before  in  the  computed  noisy 
data,  tight  clusters  occur  only  at  the  lowest  frequencies. 
The  poles  in  these  tight  clusters  are  those  which  are 
measurably  present  at  various  aspects.  The  poles  extracted 
at  higher  frequencies  are  those  which  possessed  sufficient 
measurable  energy  at  the  given  aspect.  Figure  32  depicts  the 
comparison  between  poles  extracted  from  the  measured  and 
computed  signals.  Again,  the  closest  agreement  between  the 
two  sets  of  poles  occurs  at  the  lowest  frequences. 

2.  Aircraft  Models 

Plastic  1/72  scale  aircraft  models,  coated  with 
silver,  were  used  for  transient  scattering  measurements. 
Representative  scattering  signatures  of  two  aircraft  targets, 
measured  at  six  different  aspects,  are  shown  in  Figures  33 
through  36. 

The  results  of  pole  extraction  in  target  1  are 
shown  for  a  total  of  six  different  aspects  in  Figures  37  and 
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Figure  28.  Measured  Thin  Wire  Scattering,  45  Degree  Aspect 
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Figure  29.  Measured  Thin  Wire  Scattering,  60  Degree  Aspect 
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Figure  30.  Measured  Thin  Wire  Scattering,  90  Degree  Aspect 
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Figure  33,  Target  1  Scattering,  30  Degrees  from  Nose  on 
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Figure  34.  Target  1  Scattering,  Nose  on 
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Figure  35.  Target  2  Scattering,  30  Degrees  from  Nose  on 
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38.  The  poles  extracted  at  all  six  aspects  are  shown  in 
Figure  39.  Only  one  clearly  discernible  cluster  is  present 
in  each  of  the  three  figures.  At  higher  frequencies,  no 
useful  information  is  imparted  by  the  data.  Results  of 
similar,  though  slightly  improved  quality,  were  obtained  from 
target  2.  These  results  are  presented  in  Figures  40  through 
42  in  the  format  of  Figures  37  through  39  respectively. 

Although  the  Kumaresan-Tuf ts  algorithm  is 
capable  of  extracting  low  frequency  poles  acceptably,  the 
inconsistent  results  at  higher  frequences  reveals  the  inherent 
weakness  in  an  algorithm  capable  of  processing  only  the  late¬ 
time  portion  of  a  target’s  radar  response. 

A  side-by-side  comparison  of  poles  obtained  from 
both  aircraft  by  both  the  Kumaresan-Tuf  ts  method  and  the 
Cadzow-Solomon  method  is  presented  at  the  end  of  the  chapter 
to  illustrate  the  gains  afforded  by  processing  the  early— time. 

C.  CADZOW-SOLOMON  ALGORITHM 

Recall  from  the  results  depicted  in  Figure  8  that  a  late 
transition  to  late-time,  and  the  consequent  reduction  of 
signal  power,  caused  complete  breakdown  of  the  Kumerasan-Tuf ts 
algorithm.  The  Cadzow-Solomon  algorithm  addresses  this 
shortcoming  by  processing  the  signal  at  the  instantaneous 
onset  of  early-time.  Thus,  the  Cadzow-Solomon  algorithm  is 
capable  of  processing  the  earliest  response  of  a  target  to 
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electromagnetic  excitation  where  the  response  has  the  greatest 
magnitude . 

1.  Applicability 

The  early-time  portion  of  a  target's  scattered  field 
occurs  as  long  as  there  is  a  driven  portion  of  the  total 
field.  Once  the  field  no  longer  contains  a  scattered  response 
due,  in  part,  to  the  incident  excitation  at  points  on  the 
object,  early-time  ceases  and  late— time  begins.  Hence,  the 
Cadzow-Solomon  models  both  the  system's  input  and  output,  and 
equivalently,  the  poles  and  zeros  of  the  system  transfer 
function . 

2.  Equations 

The  Cadzow-Solomon  algorithm  extends  the  auto¬ 
regressive  equation  (9),  used  in  Prony’s  method,  to  the  more 
general  autoregressive  moving  average  (ARMA)  equation 

y  =Ib|y„  +Sai^r,-.  (23) 

where  the  second  summation  term  models  the  excitation  to  the 
system . 

A  set  of  M  such  equations  in  matrix  form  is  given  by 
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As  in  the  Kumaresan-Tuf ts  method,  M  is  selected  to  be  greater 


than  the  column  dimension  of  the  data  matrix  which  is 

Kd+Kn+1  . 

3.  Excess  Poles  and  Noise  Removal 

The  Cadzow-Solomon  method  used  in  this  thesis  is  a 
modification  which  incorporates  the  non-causal  arrangement  of 
the  system  equations  used  by  Kumaresan-Tuf ts .  This 

modification  was  first  discussed  by  Norton  in  [16] .  The 
Kumaresan  approach  of  overestimating  the  system  order  can  be 
used  as  before  in  a  non-causal  model  to  constrain  the  noise 
poles  inside  the  unit  circle,  while  SVD  forces  the  signal 
poles  outside  the  unit  circle. 

Since  the  input  waveform  is  known,  its  order  can  be 
almost  exactly  determined.  In  all  the  work  of  this  thesis, 
the  input  waveform  used  is  the  double  Gaussian  depicted  in 
Figure  14.  Approximately  25  samples  defining  this  pulse  of 
0.5  nanoseconds  duration  makes  equal  25  in  equation  (23) . 
Since  the  input  is  causal,  the  signal  zeros  fall  inside  the 
unit  circle  where  they  cannot  be  easily  segregated  from 
similarly  located  noise  poles.  However,  the  signal  zeros 
impart  no  information  about  the  target  and  need  not  be 
extracted.  The  inclusion  of  the  input  in  the  data  matrix  is 
nevertheless  vital  to  the  model  of  the  system  and  the  accurate 
determination  of  the  signal  poles. 
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The  ARMA  equation  of  (23)  can  be  modified  to  obtain 


y  =5;b',y  +Sa,Xn- 
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(25) 


The  recursive  portion  of  (25)  is  now  in  a  non-causal  form 
similar  to  expression  (12)  .  A  set  of  M  such  equations  in 
matrix  form  is  given  by 
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Or,  in  matrix  notation 
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4.  Singular  Value  Decomposition 

Like  the  system  equations  of  the  Kumaresan-Tuf ts 
model,  the  system  equations  in  (26)  are  processed  using 
singular  value  decomposition.  The  coefficient  vector  is  again 
the  minimum-norm  solution,  which  constrains  the  extraneous 
poles  and  extraneous  zeros  to  be  inside  the  unit  circle. 

5.  Bias  Compensation  in  the  Cadzow-Solomon  Formulation 
By  compensating  the  eigenvalues  of  the  j;  matrix  in 

(16),  the  performance  of  the  Kumaresan-Tuf ts  algorithm  is 
significantly  improved  in  the  presence  of  noise.  Cadzow- 
Solomon  have  shown  [5]  that  if  the  actual  orders  and  K,|, 


If  the  input  and  output  noise  variances  are  not  equal,  the 
eigenvalue  shifting  theorem  used  in  Kumaresan-Tuf ts  cannot  be 
used  to  analytically  predict  the  requisite  eigenvalue 
compensation  of  Nevertheless,  when  the  input  and 
output  variances  were  assumed  equal,  and  eigenvalue 
compensation  similar  to  that  used  in  Kumaresan-Tuf ts  was 


71 


performed,  the  results  were  consistently  superior  to  those 
obtained  without  compensation.  Therefore,  the  results  of 
Cadzow-Solomon  signal  processing  presented  in  this  thesis  were 
obtained  using  eigenvalue  compensation  and  the  assumption  of 
equal  noise  variance. 

6 .  Performance 

The  Cadzow-Solomon  algorithm  was  programmed  in  Fortran 
and  tested  on  the  same  data  used  for  evaluating  the  Kumaresan- 
Tufts  algorithm.  Note  that  the  Cadzow-Solomon  algorithm 
can  use  the  early-time  portion  of  the  data  that  the  Kumaresan- 
Tufts  algorithm  can  not  use.  The  program  appears  in 
Appendix  B. 

a.  Synthetically  Generated  Data 

The  smarting  point  for  evaluating  the  performance 
of  the  Cadzow-Solomon  algorithm  was  with  synthetically 
generated  data  of  the  form  given  by  (6)  plus  the  addition  of 
input  data  required  to  model  early  time  data. 

1.  Noise  Performance 

The  algorithm  was  evaluated  at  various  signal- 
to-noise  ratios,  ranging  from  90.0  dB  to  7.0  dB .  Figure  43 
shows  the  signal  produced  by  two  s-plane  poles  at  90.0  dB , 
with  a  late-time  beginning  at  10.0  nanoseconds.  Figures  44 
through  48  depict  the  poles  extracted  from  this  signal  at  the 
different  signal-tc-noise  ratios. 


The  figures  chart  the  steady  degradation  of 
the  algorithm's  performance  with  the  increase  of  noise.  At 
30.0  dB,  the  location  of  the  low  frequency  pole  is  already 
slightly  displaced.  More  significant  is  the  location  of  one 
of  the  extracted  poles  in  the  noise  signal  space.  At  20.0  dB , 
the  low  frequency  pole  is  located  in  some  trials  on  the  real 
axis.  At  10.0  dB,  all  the  extractions  are  located  on  the  real 
axis  and  at  7.0  dB  their  locations  there  are  dispersed.  The 
extraction  of  the  higher  frequency  pole  is 
uncharacteristically  more  accurate  than  that  of  the  low 
frequency  pole.  Even  at  7.0  dB,  the  high  frequency  pole  is 
located  with  excellent  accuracy.  The  location  of  the  low 
frequency  pole  near  the  real  axis  was  chosen  deliberately  to 
illustrate  the  difficulty  in  resolving  the  slight  frequency 
difference  between  the  true  pole  and  a  noise  pole  located  on 
the  real  axis.  Also,  fewer  points  were  processed  using  the 
Cadzow-Solomon  method  than  were  processed  using  the  Kumaresan- 
Tufts  method,  since  the  largest  data  matrix  allowed  by  the 
programs  in  Appendices  A  and  B  contain  fewer  data  points  in 
the  Cadzow-Solomon  data  matrix  than  in  the  Kumaresan-Tuf ts 
data  matrix.  The  results  demonstrate  the  need  to  process  a 
substantial  number  of  points  in  order  to  accurately  extract 
low  frequency  poles. 


79 


b.  Thin  Wire  Integral  Equation  Generated  Data 

The  performance  of  the  Cadzow-Solomon  algorithm 
was  evaluated  using  the  same  set  of  data  tested  by  the 
Kumaresan-Tuf ts  algorithm.  The  results  are  presented  in 
Figure  49.  Tight  clusters  appear  at  frequencies  higher  than 
those  obtained  with  the  Kumaresan-Tuf ts  algorithm.  Figure  50 
depicts  the  poles  extracted  from  the  same  signal  at  a  20.0  dB 
SNR.  The  clustering  at  this  SNR  is  comparable  to  the  results 
obtained  by  the  Kumaresan-Tuf ts  method  with  the  noiseless 
data.  Further  angle-by-angle  comparisons  of  the  poles 
extracted  from  the  noiseless  data  and  the  20.0  dB  data  are 
depicted  in  Figures  51  through  54.  Note  the  small  number  of 
poles  in  Figure  54  due  to  the  unexcited  odd-symmetric  poles 
at  90°  aspect. 

One  further  test  was  conducted  on  computed  data 
at  a  7.0  dB  SNR.  The  results  are  depicted  in  Figure  55.  Even 
at  7.0  dB ,  discernible  clusters  are  present.  Angle-by-angle 
comparisons  of  the  poles  obtained  in  7.0  dB  data  and  those 
obtained  in  noiseless  data  are  presented  in  Figures  56  through 
59  . 

c.  Scale  Models 

The  same  scale  models  used  to  evaluate  the 
Kumaresan-Tuf ts  algorithm  were  used  to  evaluate  the  Cadzow- 
Solomon  algorithm. 


80 


Imaginary  z 


Im.igmary  z 


Extracted  poles 


Rea!  z 


Figure  51.  Integral  Equation  Thin  Wire  Comparison,  Noisel 
vs.  20.0  dB  SNR,  30  Degree  Aspect 
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Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  45  Degree  Aspect 
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Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  60  Degree  Aspect 
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Figure  54.  Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  20.0  dB  SNR,  90  Degree  Aspect 
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Figure  55.  Cadzow-Solomon  Poles,  7.0  dB  SNR 
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Figure  56. 


Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  30  Degree  Aspect 
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Figure  57.  Integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs.  7.0  dB  SNR,  45  Degree  Aspect 
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Figure  58. 


integral  Equation  Thin  Wire  Comparison,  Noiseless 
vs!  7.0  dB  SNR,  60  Degree  Aspect 
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Figure  60  depicts  the  poles  extracted  from 
measurements  of  a  0.1  meter  wire.  Three  tight  clusters  appear 
at  the  lowest  frequencies  and  at  the  highest  frequencies.  The 
noles  in  between  can  not  be  easily  discriminated.  The 
dispersion  of  these  poles  is  apparently  due  to  the  aspect 
dependence  of  their  measurable  power.  In  other  words,  these 
poles  are  excited  more  at  some  aspects  then  at  others. 

Figure  61  depicts  the  comparison  between  poles 
extracted  from  computed  data  and  measured  data.  As  in  Figure 
60,  close  agreement  exists  at  the  highest  and  lowest 
frequencies.  The  results  are  much  more  favorable  than  those 
similarly  obtained  by  the  Kumaresan-Tuf ts  algorithm. 

2.  Model  Aircraft 

Figures  62  through  64  depict  poles  extracted 
from  aircraft  target  1.  As  in  the  Kumar esan-Tuf ts  testing, 
the  Cadzow-Solomon  testing  was  conducted  at  six  different 
aspects.  Results  for  target  2  are  depicted  in  Figures  65 
through  67.  The  results  of  both  targets  show  clearly  defined 
clusters.  The  first  two  clusters  of  target  2  are 
exceptionally  tight.  However,  the  mid-frequency  clusters  of 
target  2  are  not  as  clearly  formed  as  those  of  target  1. 

Comparisons  of  poles  obtained  with  each  method 
for  target  1  and  2  are  depicted  in  Figure  68  and  69 
respectively.  These  two  figures  graphically  depict  the  clear 
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Figure  60. 


Cadzow-Solomon  Poles,  Measured  Thin  Wire 
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Figure  61.  Thin  Wire  Comparison 
Equation 
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Figure  62.  Cadzow-Solomon  Poles  Target  1,  Three  Aspects 
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1  Figure  63.  Cadzow-Solomon  Poles  Target  1,  Three  Aspects 
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Cadzow-Solomon  Poles  Target  2,  Three  Aspects 
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Figure  69.  Pole  Comparisons,  Target  2,  All  Six  Aspects 


superiority  of  the  Cadzow-Solomon  algorithm  over  the 
Kumaresan-Tuf ts  algorithm. 

In  order  to  obtain  an  initial  indication  of 
the  possibility  for  target  classification  through  pole 
extraction,  nose-on  measurements  of  two  additional  aircraft 
models  were  made,  processed  and  compared  with  the  results  of 
targets  1  and  2.  The  nose-on  measurements  of  targets  3  and 
4  appear  in  Figures  70  and  71  respectively.  A  comparison  plot 
of  poles  extracted  from  each  of  the  four  targets  is  depicted 
in  Figure  72.  Each  of  the  four  aircraft  measured  are  fighters 
of  similar  size  and  shape  (see  Table  1).  The  poles  for  each 
target  are  sufficiently  different  in  this  single  measurement 
to  identify  each  aircraft  individually.  However,  some  of  the 
poles  are  arranged  in  clusters  which  appear  with  a  harmonic 
pattern  similar  to  that  obtained  for  either  of  the  first  two 
aircraft  at  various  aspects.  In  order  to  more  fully  assess 
the  target  classification  capability  of  pole  extraction, 
several  measurements  should  be  made  of  a  given  aircraft  model. 
A  plot  of  the  poles  extracted  from  each  of  these  measurements 
would  form  clusters  at  the  locations  of  the  true  poles.  The 
centroid  of  each  of  these  clusters  would  then  be  compared 
against  the  centroid  poles  similarly  obtained  from  other 
aircraft.  Although  several  poles  of  different  aircraft  might 
be  similar,  the  set  of  poles  belonging  to  an  aircraft  could 
form  the  basis  for  classification  if  that  set  was  unique  among 
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Nose  on 


Time  (nanoseconds) 


Figure  70.  Target  3  Scattering,  Nose-on 


104 


Nose-on 


Figure  71.  Target  4  Scattering,  Nose-on 
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the  sets  belonging  to  all  other  measured  aircraft.  The  results 
in  Figure  72  demonstrate  the  possibility  of  using  the  Cadzow- 
Solomon  pole  extraction  algorithm  to  aid  in  the  classification 
of  aircraft,  perhaps  by  use  of  the  extracted  poles  in 
constructing  annihilation  filters. 

TABLE  1.  FULL  SIZE  DIMENSIONS  OF  TARGETS  RECORDED 


Target  number 

1 

2 

3 

4 

Overall  length 

12 . 20 

15.03 

16 . 94 

16.00 

(meters ) 

Overall  height 

3 .35 

5.09 

4 . 51 

4 .80 

(meters ) 

Wingspan 

10.96 

10.00 

11 . 43 

13.95 

(meters ) 


Tailplane  span 
(meters ) 


Unknown 


5.58 


6.92 


5.75 


III.  SUMMARIES  AND  CONCLUSIONS 


In  this  chapter,  a  step-by-step  guide  through  each 
algorithm  is  presented.  At  each  step,  techniques  and  lessons 
learned  are  discussed  together  with  general  observations. 
Conclusions  are  presented  at  the  end  of  the  chapter. 

A.  KUMARESAN-TUFTS 

The  first  step  in  processing  a  signal  with  the  Kumaresan- 
Tufts  algorithm  is  to  determine  the  beginning  of  early-time. 
The  objective  is  to  pick  the  earliest  possible  starting  point 
without  entering  into  the  latter  part  of  early-time.  If  the 
starting  point  for  processing  is  improperly  chosen  to  include 
the  early-time,  the  results  will  be  completely  unreliable 
since  the  signal  no  longer  satisfies  the  late  time  model.  If 
the  starting  point  is  chosen  too  late,  the  signal  may  not  be 
sufficiently  strong  in  the  presence  of  measurement  noise. 
Since  the  signal  is  the  sum  of  exponentially  damped  sinusoids, 
the  optimum  starting  point  is  at  the  precise  instant  of 
transaction  into  late-time.  The  key  to  determining  the 
beginning  of  late-time  is  in  determining  the  beginning  of 
early  time.  Determining  the  first  response  of  the  target  to 
excitation  cannot  usually  be  done  by  a  simple  visual 
inspection  of  measurement  data.  Unless  the  exact  distance  to 
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the  target  is  known,  the  most  accurate  method  attempted  by  the 
author  for  determining  the  beginning  of  early— time  is  to 
process  the  signal  using  the  Cadzow-Solomon  algorithm.  This 


is  discussed  in  the  next  section.  However,  the  reliance  of 
the  Kumaresan-Tuf ts  algorithm  on  information  provided  by  the 
Cadzow-Solomon  algorithm  is  an  obvious  disadvantage  of  the 
former  method. 

Once  the  starting  point  for  processing  has  been  selected, 
the  next  step  is  to  determine  the  dimensions  of  the  data 
matrix  and,  consequently,  the  number  of  points  in  the  signal 
to  be  processed.  In  trials  conducted  on  noiseless  synthetic 
data,  the  accuracy  of  pole  extraction  increased  steadily  with 
the  increase  in  the  data  matrix  dimensions.  These  trials  were 
conducted  up  to  the  limit  of  the  array  dimensions  defined  in 
the  computer  program  of  Appendix  A.  The  number  of  points 
processed  in  measurement  data  should  be  as  large  as  possible, 
while  still  meeting  the  following  two  constraints.  First, 
incorporate  as  many  cycles  of  the  data  as  possible.  Usually, 
visual  inspection  of  the  data  reveals  a  repeating  pattern 
which  should  be  entirely  incorporated  into  the  window  of 
points  to  be  processed.  When  only  portions  of  these  patterns 
are  selected,  a  disproportionate  weighting  tends  to  be  placed 
on  certain  poles.  Second,  signal  portions  late  in  the 
response  which  are  no  longer  distinguishable  in  the  presence 
of  measurement  noise  should  not  be  selected. 
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The  final  step  involves  determining  the  number  of  true 
poles  in  the  system.  The  following  approach  has  proven  to  be 
the  most  successful.  First,  process  the  signal  without  any 
eigenvalue  compensation  to  establish  an  upper  bound  on  the 
order  of  the  system.  In  most  cases,  the  number  of  poles 
outside  the  unit  circle  will  be  less  than  the  overestimated 
order  of  the  system.  If  increase  the  row  dimension  of 
the  data  matrix  in  order  to  increase  the  estimated  order  of 
the  system,  and  repeat.  When  the  I'umber  of  poles  is  less  than 
the  estimated  order  of  the  system,  then  one  should  gradually 
increase  the  number  of  eigenvalues  compensated  in  successive 
trials,  while  closely  observing  the  effects  induced  on  the 
poles  outside  the  unit  circle.  As  the  number  of  eigenvalues 
compensated  is  steadily  increased,  noise  poles  and  weak  signal 
poles  will  move  inside  the  unit  circle.  The  programs  in 
Appendix  A  and  B  allow  the  user  to  compare  the  results  of 
successive  trials,  by  generating  overlays  for  each  plot.  If 
N  poles  are  in  the  signal  space,  at  least  the  first  N 
eigenvalues  must  not  be  compensated,  or  true  poles  may  be 
lost.  As  the  actual  order  of  the  system  is  approached  by 
compensation,  the  user  will  notice  an  orderly,  even 
arrangement  assumed  by  the  noise  poles.  If  certain  poles 
still  remain  suspect  after  compensation,  vary  slightly  the 
other  parameters,  such  as  the  starting  point  and  the 
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dimensions  of  the  data  matrix.  Generally,  only  true  signal 
poles  will  repeatedly  assert  themselves  under  varying 
parameters . 

B.  CADZOW-SOLOMCN 

The  techniques  and  general  observations  offered  in  the 
preceding  section  apply  equally  to  the  Cadzow-Solomon 
algorithm.  An  important  consideration  in  this  method,  not 
discussed  above,  is  the  selection  of  the  beginning  of  early- 
time.  Candidates  for  a  starting  point  are  usually  at  or  near 
zero  crossings  within  approximately  thirty  points  of  the 
object’s  first  definite  response  to  electromagnetic 
excitation.  Begin  processing  ac  the  chosen  point  while 
varying  parameters  in  successive  trials.  Select  the  point 
whose  successive  results  are  the  most  consistent  under  varying 
parameters . 

The  selection  of  the  starting  point  for  beginning  of 
early-time  can  be  very  critical.  For  example,  not  a  single 
pole  could  be  extracted  in  one  trial  wherein  the  starting 
point  occurred  only  ten  points  after  the  actual  starting 
point.  Additionally,  in  most  cases  observed,  the  late-time 
start  given  by  the  selected  early-time  occurred  within  less 
than  two  points  from  a  zero  crossing.  If  this  observation 
proves  to  be  generally  true  in  later  research,  it  may  serve 


1 
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as  a  way  to  check  the  starting  point  selected  for  one 
algorithm  in  terms  of  the  other. 

C.  CONCLUSIONS 

Both  the  Kumaresan-Tuf ts  and  the  Cadzow-Solomon  algorithms 
can  effectively  extract  poles  from  the  scattering  response  of 
a  radar  target.  Because  both  algorithms  obtain  a  least- 
squares  solution  to  the  system  model,  both  perform  acceptably 
in  the  presence  of  noise.  Although  eigenvalue  compensation 
is  not  analytically  justified  in  the  Cadzow-Solomon  algorithm, 
the  results  obtained  through  eigenvalue  compensation  in  this 
method  were  generally  superior  to  those  similarly  obtained  in 
the  Kumaresan-Tuf ts  method.  The  results  demonstrated  the 
inherent  advantages  of  an  algorithm  capable  of  processing  a 
target's  strongest  response  in  the  early  time.  The  Kumaresan- 
Tufts  method  compared  favorably  with  the  Cadzow-Solomon  only 
in  responses  with  a  long  late— time. 
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APPENDIX  A. 


THE  KUMARESAN-TUFTS  POLE  EXTRACTION  ALGORTITHM 


The  following  program  implements  the  Kumaresan-Tuf ts 

algorithm  as  described  in  Chapter  2  of  this  thesis.  The 

program  is  written  in  Fortran  77.  The  SVD  and  root-finding 

subroutines  called  by  this  program  are  found  in  the  EISPACK 

library  [18] .  The  SVD  subroutine  is  a  translation  from  ALGOL 

as  given  in  [19]  .  The  matrix  multiplication  and  graphics 

subroutines,  also  called  by  this  program,  are  found  in 

Appendix  C  and  D  respectively. 

INTBCaR  n»R,Kd,M,)«,MAGKL,NSTWPT,iaLTAY 
INmaF  m,NCAUS,l«ENU,L/l/ 

INra3R*2  KdPLT 

RERL*8  A(70,70),W(70),U(70,70),V{70,70),RV1(70) 

RERL*8  VS (70, 70) ,OT{70,70) ,AINV{70,70) ,X(70) 

REAL*8  XP(70)  ,B{70)  ,Siaffi(70,70)  ,SIG(70,70) 

REAL*8  aF(70),RO0TR(70),ROOTK70) 

REAL*8  D(1024),AVG,MAanP/1.0E-16/,Dy(140) 

0CHPIiX»16  S(70) 

LOGICAL  MATO/. TRUE./, MATV/.'mUE./,CAUSAL/.TRUE./,LCIIG/.TOJE./ 

LOGICAL  DSEr/.FALSE./,NUFILE/.TRUE./ 

CHARACTER  TrnJC*16,HEAISR*64,YN*l,DC*l,TmJER*16,TrnZI*16 
CHARACTER  TnL*16 


Eliter  parameters  for  processing 


11  IF  {KET)  CLOSE  (10) 
M0V!RLAY=O 
OPEN  (10,  FIIi)=’ PLOT’) 
IF  (DSET)  00  TO  85 


WRITE 
WRITE 
WRITE 
WRITE 
WRITE  (* 
WRITE  {* 
WRITE  (* 
WRITE  (* 


Vfelootie  to  signal  processing  using  the' 
Kunaresan-Tiifts  method’ 

I 

Do  you  want  ' 

1 

1.  The  long  version  for  beginners' 

2.  The  short  version  for  pros' 
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15  WRITC  (*,*)  'Please  enter  1  or  2  ’ 
READ  {*,*)  N 
IF  (N  .m.  1)  ■IBEN 
LCN&=,‘reUE. 

lUEIF  (N  .EQ.  2)  'nON 
LCM&=.FAliSE. 

BjSE 
GO  TO  15 
EMDIF 


16 


10 


WRITE  (*,*)  'Session  will  begin  with  entry  of  paraneters  needed  fo+r  processing' 
WRITE  (*,*) 


WRITE  (*,*)  'Do  you  want  to  enter  parameters  fran' 
WRITE  (*,*)  '  ' 


WRITE  (*,*)  '1.  The  te^xard' 

WRITE  {*,*)  '2.  A  previously  created  file  of  parameters' 
WRITE  (*,*)  '  ' 


WRITE  (*,*)  'Please  enter  1 
ROD  (*,*)  N 
IF  (N  .03.  1)  TWW 


or  2 


GO  TO  1 

ILSEIF  (N  .03.  2)  TWIf 

WPITE  (*,*)  'Biter  title  of  file  containing  parameters' 

READ  (*,105)  TTIL 

0PBJ{l,FIIi5=Tm) 

READ(1,105)  Tmi: 

READ (1,110)  NPTS 
READ(1,110)  MfT 
READ (1,110)  Kd 
READ (1,110)  M 
READ  (1,110)  KLTAY 
READ (1,110)  NSTETTT 
READ(1,110)  NCAUS 
CLOSE(l) 

00  TO  85 
ILSE 
00  TO  16 


□IDIF 

WRITE  (*,*)  '  ' 


1  NUFII£=.TTOE. 

IF  (.HOT.  DSET)  NSTRTPT=1 

WRITE  (*,*)  'Biter  title  of  data  file  to  be  read' 

READ  (*,105)  TTTTE 
CPDJ(l,FIIiE=TTn£) 

READ(1,105)  HOUSE 
READ (1,110)  NPTS 
IF  (NPTS  .or.  1024)  THTO 

WRITE  (*,*)  '(Amber  of  points  in  data  file  eso^eds  the  dimension' 
WRITE  (*,*)  'of  the  array  used  in  the  program  to  store  the  file' 
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STS' 

H©IF 

CL06E(1) 

IF  (DSET)  TBEN 

IF  OBreTPTf(Kd+M-l)*ISIiTRY  .l£.  NPTS)  GO  TO  85 
EXDIF 

3  IF  (NUFIIi:)  TOEJJ 

Wirre  (*,*)  'Biter  Kd,  >=  the  estimated  order  of  the  systan  ' 
REAP  (*,*)  Kd 
IF  (Kd  .GT.  69)  TOBJ 

HRTIE  (*,*)  'Kd  must  be  less  than  70,  or  dimensicn  statements' 
WRITE  (*,*)  'in  this  program  must  changed  by  the  user' 

00  TD  3 

ILSEIF  (Kd  .LT.  2)  TTOl 

WRITE  (*,*)  'Kd  must  be  at  least  2' 

GO  TO  3 
HffilF 

IF  {2*Kd  ,Gr.  NPTS)  TaOi 

WRITE  {*,*)  'Kd  must  be  less  than  or  equal  to  ',NPTS/2 
GO  TO  3 

E'.SZIF  (2*Kd  ,B2.  NPTS)  'lUEN 
WRITE  (*,*)  'Kd  equals’, Kd 
WRITE  (*,*)  'M  must  be',Kd 
WN(d 

WRITE  (*,*)  'since  there  are  a  total  of, NPTS 

WRITE  (*,*)  'points  in  '.TITUE 

00  TO  45 

FNDIF 

GO  TO  4 

ELSEIF  (DSET)  THEN 
N=M 

20  IF  (NSTRTPT+(N+M-l)*raLTAY  .L£.  NPTS)  THBi 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

25  WRITE  {*,*)  'Kd  may  range  from  ’,NRT 

MOTE  (*,*)  '  to',N 

WRITE  (*,*)  'Biter  Kd' 

READ  (*,*)  Kd 

IF  (Kd  .<31.  NRT  .AND,  Kd  .!£.  N)  00  TO  85 
00  TO  25 
ELSE 
N=N-1 
00  TO  20 
IMDIF 
EUDIF 

4  IF  (NUFII£)  THEN 

WRITE  (*,*)  'Ehter  M,  the  row  dimension  of  the  data  matrix' 

IF  (.NOT.  DSET  .AND.  LONG)  THEN 
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MRITC  (*,*)  '  ' 

WRITC  (*,*)  'Note:  Kfl+M  points  in  title 
WRITE  (*,*)  '  will  be  processed  ' 

WRITE  (*,*)  ’  ' 

ENDIF 

30  WRITE!*,*)  'M  may  range  from',  Kd 
IF  (NPTS-Kd  .GT.  69)  THfK 
WRITE  (*,*)  '  to  69' 

BjSE 

^’RTIE  (-,*)  '  to', NPTS-Kd 

QIDIF 

READ  (*,*)  M 
IF  (M  .CT.  69)  TTIEM 

WRITE  (*,*)  'M  must  also  be  less  than  70* 

00  TO  30 

FLSEIF  {M  .LT.  Kd)  THQ! 

WRITE  (*,*)  'M  must  be  greater  than  or  equal  to  Kd,  Kd=  ',Kd 
GO  TO  30 

OiSEIF  (Kd-W  -GT.  NPTS)  THEN 

WRITE  (*,*)  'Kd+M  must  be  less  than  or  equal  to',NPrs,',’ 
WRITE  (*,*)  'the  number  of  data  points  in',Trni: 

WRITE  (*,*)  '  ' 

GO  TO  30 
ENDIF 
EIjSE 
IHCd 

35  IF  (I*STTmTV(Kd-W-l)*rCLTAY  .LE.  NPTS)  TEEN 
N=Nfl 
GO  TO  35 
ELSE 
N=N-1 
ENDIF 

IF  (N  .EC.  Kd)  TEEN 

WRITE  (*,*)  'Mmust  equal’, Kd 

EHCd 

GO  TO  85 

ENDIF 

IF  (N  .GT.  69)  N=69 

40  WRITE  (*,*)  'M  may  range  from',Kd 

WRITE  (*,*)  '  to',N 

WRITE  {*,*)  'Ehter  M' 

RIM  (*,*)  M 

IF  (M  .($.  Kd  .AND.  M  .!£.  N)  GO  TO  85 

00  TO  40 

ENDIF 

45  IF  (.NOT.  NUFI1£)  00  TO  85 

5  N=1 

50  IF  (NSTETPT+N*(Kd-Hi-l)  .1£.  NPTS)  TEEN 
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KWH-1 
GO  TO  50 
OiSE 
N=W-1 
EHDIF 

IF  (N  .EQ.  1)  IHEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WtlTE  (*,*)  'Spacino  can  only  be  1' 

ismY=i 

IF  OWFIUl)  THEN 
GO  TO  50 
ELSE 
GO  TO  85 
E>©IF 
EM)IF 

IF  (.NOT.  DSET  .W©.  LONG)  THEN 

WRITE  (*,*)  'Ehter  spacing  between  the  ',K(Wf 

WRITE  (*,*)  'data  points  of  '.TTTTZ 

WRITE  (*,*)  'to  be  processed  ' 

WRITE  (*,*)  '  ’ 

WRITE  {*,*)  'If,  for  exanple,  one  is  chosen,  then  ',Kd+M 
WRITE  (*,*)  'consecutive  points  in  ',TITT£ 

WRITE  (*,*)  'will  be  processed  ' 

WRITE  (*,*)  '  ' 

EiroiF 

55  WRITE  {*,*)  'Spacing  may  range  frcn  1  ' 

WRITE  (*,*)  '  to’,N 

READ  (*,*)  EELTAY 

IF  (DELTAY  .OT.  1  .AND.  OELTAY  .LE.  N)  TEEN 

IF  (nuehz)  then 
GO  TO  60 
ELSE 
GO  TO  85 
ENDIF 
ELSE 
GO  TO  55 
EKDIF 

60  WRITE  (*,*)  'Do  you  wish  to  adjust  eigenvalues?  (y/n)’ 

READ  (*,120)  YN 

IF  (YN  .E)3.  'N'  .OR.  YN  .EJ2.  'n')  THEN 
IF  (nufuje)  go  To  6 
GO  TO  85 
ENDIF 

IF  (YN  ,NE.  'Y'  .AND.  YN  .NE.  'y')  GO  TO  60 
2  WRITE  (*,*)  'Discard  or  oonpensate  eigenvalues?  (d/c)' 

READ  (*,120)  DC 

IF  (DC  .EC.  'D'  .OR.  DC  .EC-  'd')  GO  TO  65 
IF  (DC  .NE.  'C  .AND.  DC  .NE.  'c')  GO  TO  2 
WRITE  (*,*)  'Ehter  estimate  of  the  actual  order  of  the  system' 
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WRITE  (*,*)  '  ' 

IF  (liCNG)  THEN 

WRITC  (*,*)  'This  estimate  will  be  used  to  detenniae  the  ' 
WRITE  (*,*)  'mmber  of  eigenvalues  conpensated  or  discarded  ' 
BIDIF 

65  WRITE  (*,*)  'the  estimate  may  range  froB  2' 

WRITE  (*,*)  '  to',Kd-l 

REM)  (*,*)  NRT 

IF  (WT  .or.  Rd  .OR.  Mrr  .LT.  2)  TT©) 

GOTO  65 

EliSEIF  (.NCrr.  NUFn£)  THEN 

00  TO  85 

QIDIF 


6  MSTHTPT=1 

70  IF  (N?IHTPT+(Kd+M-l)*r*LTAY  .!£.  NPK)  THEN 

NSTim>T=WSTHTmi 
GO  TO  70 
ELSE 

NSTHTT>T=«STHTPT'l 

HJDIF 

IF  (NSTHTPT  .E)3.  1)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 
WRITE  (*,*)  'the  starting  point  for  processing  the  data' 
WRITE  (*,*)  'must  be  the  first  point  in  the  data  file' 

00  TO  85 
ENDIF 

WRITE  (*,*)  'Enter  desired  starting  point  in  data  file’ 

IF  (.NOT.  DSET  .AND.  I£NG)  THEN 

WRITE  (*,*)  '1  indicates  the  first  point  in  the  data  file  ' 
ENDIF 

WRITE  (*,*)  '  ' 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far, ' 

75  WRITE  (*,*)  'the  starting  point  may  range  from  1' 

WRITE  (*,*)  '  to',NSTHTHT 

REM  (»,*)  N 

IF  (N  .GE.  1  .AND.  N  .!£.  NSTRTHT)  THEN 

KSTHTPT=N 

ELSE 

WRITE  (*,*)  'Ehter  starting  point  again’ 

WRITE  (*,*)  '  ' 

00  TO  75 
ENDIF 

IF  (.NOT.  NUFI1£)  00  TO  85 
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WRITE  (*,*)  'Do  you  want  the  data  matrix  arrangement  to  be' 
WRITE  (*,*)  '  ' 

WRITE  (*,*)  '1.  Causal' 

WRITE  (*,*)  '2.  Non-causal’ 

WRITE  (*,*)  '  ' 
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80 


WRITC  (*,*)  'Please  alter  1  or  2  ' 

READ  (*,*)  MCAUS 
IF  (HCMJS  .US.  1)  THEU 
C3U]SAIi=.'naJE. 

naeiF  oous  .bq.  2)  iben 

aUSAL=.FALSE. 

ELSE 
GO  TO  80 
ENDIF 
00  TO  85 

9  KRTTZ  (*,*)  'filter  title  of  file  to  contain  parameters' 
READ  (*,105)  TITL 
Orai(l,FTI£=TnL) 

WRTIEd.lOS)  TmE 
WRITE (1,110)  NPTS 
WRITE (1,110)  NRT 
WRITE  (1,110)  Kd 
WRITE  (1,110)  M 
WRnE(l,110)  KLTAY 
WRITE  (1,110)  NSTRTPT 
WRnE(l,110)  NCAUS 
CLOSE(l) 

IF  (DSET)  00  TO  85 

12  IF  (DSET)  THEN 
CL0SE(2) 
diOSEO) 

CALL  SUBPLT(NDVERIiAY) 

EJIDIF 


85 


DSin'=.TRUE. 


NUFI1£=.FALSE. 

WRnE(*,*) 

'  ' 

WRnE(*,*) 

'1. 

+rni: 

WRnE(*,*) 

f 

WRnE(*,*) 

'2. 

WRnE(*,*) 

'3. 

WRITE(*,*) 

'4. 

WRnE(*,*) 

'5. 

+Y 

WRnE(*,*) 

'6. 

4PT 

WRnE(*,*) 

1 

-tPTHKd+E!-! 

IF  (NCAUS  , 

.E>3. 

WRnE(*,*)  '7. 
■KBAL 
ELSE 


Data  file  to  be  processed  ' ,T 

Nunber  of  data  points  in  data  file  ' ,NPTS 

Estimated  order  of  the  system  ' ,NRT 

Kd,  the  nunber  of  ooluntis  in  the  data  matrix', Kd 
H,  the  nunber  of  rows  in  the  data  matrix', M 
Spacing  between  data  points  being  processed  ' ,  DELTA 

First  point  in  the  data  file  to  be  processed ',NS1RT 

Last  point  in  the  data  file  to  be  processed ',NSTRT 


1)  TEEM 

Data  matrix  arrangenent  for  processing  CA 


IF 
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WRrn:(*,*)  n.  Data  matrix  arrangeoent  for  wocessing  NCN-CA 
■KJSAL  ' 

EXDIF 

’  ■ 

WRnT;(*,*)  '8.  Begin  processing  using  above  settings' 

WOTE:(*,*)  '9.  Store  paraneters  1-7  in  a  file’ 

WRITE(*,*)  '10.  Retrieve  parameters  1-7  frcm  a  previously  created 
+fiie’ 

WrnE(*,*)  'll.  Reset  overlays' 

WnE(*,*)  '12.  Re-plot  overlays' 

WOTEC*,*)  '13.  Bid  this  session  of  lOmaresan-Tufts  signal  process 
+ing' 

VRITE{*,*)  ’  ' 

MRrnK*,*)  'Biter  an  integer  from  1  to  12  to  make  dianges  as  often 
+  as  you  desire' 

90  READ  (*,*)  »1B«J 

IF  (»®RJ  .LT.  1  .OR.  MffiMU  .CT.  13)  17©} 

WRrrc(*,*)  'Biter  an  integer  from  1  to  13' 

00  Tt)  90 
fUDIF 

GO  TO  (l,2,3,4,5,6,7,8,9,10,ll,12,13),NKa}U 

8  aPBJd.FIliOTIU:) 

RERD{1,105)  HEWJTO 
READ (1,110)  NPTS 
RlSyO  (1,115)  X2 
READ(1,115)  XQ 
DO  95  1=1, NPTS 
READ(1,115)  D(I) 

95  OCWnNUE 

dJOSEd) 

KdPL'MOd 

WRTIIK*,*)  'Biter  title  of  file  to  contain  real  part  of  poles' 

RnD(*,105)  TTIUR 

C*©((2,file=TnTiR) 

VRrn:(*,*)  'Biter  title  of  file  to  contain  imaginary  part  of  poles' 
READ(*,105)  TTlUa 
CPQ}(3,file=nTl£I) 
wRrn:do,ioo)  (rjplt) 
wrre(io,i05)  TnuR 
wRrn:(io,io5)  Tnm 
100  raWAT(I2) 

»HttX(M,Kd) 

105  roRMAT(A) 
no  roRMAT(I5) 

115  raRMAT(E12.6) 
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I 

> 


120  P3RMM'(A1) 

C  Form  data  matrix 

DO  126  I=l,Kfl+« 

Dy(I)=®  ( (I-l) ‘DBL.’mY+NSTOTTT) 

125  OKnNUE 

130  DO  140  1=1, M 
DO  135  J=l,Kd 
A(I,J)=Cir(I-KJ) 

135  OOfTINUE 
140  OCNTINUE 

B(l)=Dy{l) 

DO  145  1=2, H 
B(I)=A(I-1,1) 

145  OQWriNUE 

C  Begin  singular  value  deocnqposition 

CAa  SVD{MACHn',M,Kd,MN,A,W,MATO,U,lftTV,V,im,RVl) 

C  Errors  in  SVD? 

IF  (lERR  .GT.  0.0)  TOEN 

WRITE  (*,*)  'Error  in  singular  value  number  •,in^R,STO 
HffilF 

IF  (YN  .B3.  'N')  GO  TO  190 

DO  150  1=1, Kd 
XP(I)=0.0 
150  OONTINUE 

C  Discard  or  conpensate  eigenvalues 

C  Order  singular  values 

XP(1)=M/(1) 

DO  165  1=2, Kd 
DO  160  «J=1,I 

IF  (W(I)  .GT.  XP{J))  TTOJr 
DO  155  K=I+1,J,-1 
155  XP{K)=XP(K-1) 

XP(J)^(I) 

00  TO  165 
ENDIF 

160  OCNTINUE 

XP(I+1)=W(I) 

165  OOOTINUE 

C  XP(  )  now  contains  ordered  singular  values-XP{l)  is  the  largest 
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C  Discard  eigenvalues 

IF  (DC  .!C.  'DM  TffiH 
DO  170  J^ffiTfl.Kd 
170  W(J)  =  (0.0) 

ELSE 

C  Conpensate  eigenvalues 

kVG=O.0 

DO  175  J=*«T+l,Kd 
AVG=AVG+XP(J)**2 
175  CONTINUE 

IF  (Kd  -GT.  MW)  AVG=AVG/DBIi:(FIiaAT(Kd-MW)) 

DO  185  0=1, Kd 
DO  180  K=l,Kd 

IF  (  W(J)  .EQ.  XP(K)  )  THEN 
IF  (  K  .CT.  NRT  )  THEN 
W(J)=0.0 
ELSE 

W{J)=DSi2RT{nABS(  W(J)*W{J)-AW)) 

ENDIF 
GO  TO  185 
ENDIF 
180  CONTINUE 
185  CONTINUE 
ENDIF 

190  DO  200  1=1, M 
DO  195  J=1,M 
OT(I,J)  =  (U(J,I)) 

195  CONTINUE 
200  CONTINUE 

c  Fcam  SIGMA+  (KbM) 

DO  210  1=1, Kd 
DO  205  J=1,M 
SIGMA(I,J)=0.0 

IF  (I  .EC.  J  .AND.  W{J)  .NE.  0.0)  THEN 
SIGMA(I,J)=1.0D0/W{J) 

ELSE 

SIGMA(I,J)=O.OdO 

ENDIF 

205  CONTINUE 
210  CCKTINUE 

C  Tam  SIGMA  (IM(d) 

DO  220  1=1, M 
DO  215  J=l,Kd 
SIG(I,J)=0.0 

IF  (I  .EC.  J)  SIG(I,J)=#(J) 
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215  OCWTINUE 
220  OCWriHUE 


C  \NUxKd,SIGMA-H=«dxM,VS=»axM 

caii,  iDatJL(v,si(aft,Kd,Kd,M,vs) 

C  VS=*M,lTMhM,AIN\Na^ 

OIL  M3MJL(VS,m',Ka,M,M,AINV) 

C  Calculate  matrix  multiplication  of  AINV  x  B,  where 

C  AINVNCdxM,B=«3d,XP=Kdxl 

CALL  MXMJL{AINV,B,Kd,M,L,XP) 

C  Calculate  autoregressive  coefficients  fran  prediction  coefficients 

IF  (XP(Kd)  .BS.  0.0)  THEH 
WRITE  (*,*)  'IiyJOR,  avoiding  division  by  zero' 

STO> 

QjSE 

B(Kd)=1.0d0/XP{Kd) 

EUDIF 

DO  225  1=2, Kd 
B{I-l):^(Kd)*XP(Kd-I+l) 

225  OOtmUUE 

DO  230  1=1, Kd 
X(I)=-B{Kd-I+l) 

IF  (NCAUS  .02.  1)  X(I)=-XP(Kd-I+l) 

230  OCWTINUE 
X(Kd+l)=1.0 

C  Confute  the  roots  of  the  polynctnial  in  z 

CAIi,  PCXin'(X,Ct:»',KD,ROCTO,ROCTI,IER) 

IF  (IIF  .NE.  0)  WRITE  (*,*)  'EFROR  with  PCIFT,  IIF=' ,nF,STM> 

DO  235  1=1, Kd 
WRITE(2,115)  ROCriR(I) 

WRITE(3,115)  ROOn(I) 

S (I)  =DC}lPU(ROCfrR (I)  ,ROOn (I) ) 

235  CCWriNUE 

MAGPCXrO 
DO  240  1=1, Kd 

IF  (CI»BS{S(I))  .GE.  l.OdO)  MAGPtXi=«AGPCL+l 

240  ocwnuuE 

WRrn:(*,*)  '|  of  poles  with  magnitude  <=  l',Kd-«AGPCL 
WRITE  (*,*)  'HIT  ANY  KEY  ■TO  OCWTINUE’ 

READ  (*,105)  HEADER 


1?3 


C  Plot  poles 


NC^mAY=«CVISLAY+l 

CL0SE(2) 

CL0SE(3) 

Caii  SUBPLT(NCfimAY) 

J=0 

K=0 

DO  245  1=1, Kd 

IF  {CI»BS(S(I))  .LT.  1.0)  TODJ 

J=J+1 

K=«+l 

WRITE  (*,*)  S(I),CDABS(S(I)) 

IJIDIF 

IF  (J  .5J.  20)  THH< 

WRrre  (*,*)  'Eiiter  any  )cey  to  continue* 

RIM  (*,105)  HIMIK 

J=0 

UOIF 

245  OCWriNUE 

WRrnE(*,*)  'Poles  with  magnitude  less  than  one:  ',K 

GO  TO  85 
13  m> 

ODD 
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APPENDIX  B:  THE  CADZOW- SOLOMON  POLE  EXTRACTION  ALGORITHM 


The  following  program  implements  the  Cadzow-Solomon 
algorithm  as  described  in  Chapter  2  of  this  thesis.  The 
program  is  written  in  Fortran  77.  The  SVD  and  root-finding 
subroutines  called  by  this  program  are  found  in  the  EISPACK 
library  [18] .  The  SVD  subroutine  is  a  translation  from  ALGOL 
as  given  in  [19].  The  matrix  multiplication  and  graphics 
subroutines,  also  called  by  this  program,  are  found  in 
Appendix  C  and  D  respectively. 


SLARGE 

INTTX3R  IIRR,Kd,Kn,M,MN,MA(ya,,NSTRTPT,I«LTRY 
INTBt3R  IIR,NCAUS,AMENL',INSlTm>T 
IIfrB®R*2  KdPLT 

REAL*8  A(7:,70) ,W(70) ,U(70,70) ,V(70,70) ,RV1{70) 

REAL*8  VS  (70, 70)  ,IT(70,70)  ,AINV(70,70)  ,X(70) 

REAL*8  XTC^O)  ,B(70)  ,SiaiA(70,70)  ,SIG(70,70) 

REAL*8  aF(70),ROCTO(70),RO0TI(70) 

REAL  MAG 

REAL*3  D(1024) ,AVG,MACHlP/1.0E-16/,Dy(140) ,nx(1024) 

0CMPIiX*16  S(70) 

LOGICAL  MATO/ .TRUE. / , MATV/ .TRUE. /  .CAUSIAL/ .TRUE. / , LONG/ .TRUE. / 
LOGICAL  DSET/. FALSE. /.MJFIIi:/. TRUE./ 

CHARACTER  TnLE*16,HEADER*64,YN*l,DC*l,TnLER*16,Tni£I*16 
cHARAcnR  'Tivie.TmD^ie 

C  filter  parameters  for  processing 

14  IF  (DSIT)  CLOSE (10) 

NCViRLAY=0 

CP0J(lO,niJE>=’PLOT') 

IF  (DSET)  00  TO  215 

WRITE  (*,*)  'Welcone  to  signal  processing  using  the' 

WRITE  (*,*)  'Cadzcw-Solonon  method' 

WRITE  (*,*)  ■  ' 

WRITE  (*,*)  'Do  you  want  ’ 


WRITE  (*,*) 


WRITE  (* 
WRITE  (* 
wRrrc  (* 
WRITE  (* 


1.  “nie  long  version  for  beginners' 

2.  Tbe  short  version  for  pros' 

I 

Please  enter  1  or  2  ' 


READ  (*,*)  N 
IF  (N  .BQ.  1)  TBEH 
I£N5=.reUE. 

EZiSEIF  (N  .EQ.  2)  TROT 
Lai&=.FALSE. 

OiSE 
GO  'ro  25 
HOIF 


WRITE  (*,*) 
+r  processing 
WRITE  (*,*) 


Session  will  begin  with  entry  of  parameters  needed  fo 


WRITE  {*,*) 
WRITE  (*,*) 
WRITE  (*,*) 
WRITE  (*,*) 
WRITE  (*,*) 
WRITE  (*,*) 


Do  you  want  to  alter  parameters  from' 

( 

1.  The  keyboard' 

2.  A  previously  created  file  of  parameters' 
Please  enter  1  or  2  ' 


READ  (*,*)  N 
IF  (N  .02.  1)  TTEJ 
00  TO  8 

ELsm  (N  ,EQ.  2}  mu 

WRITE  (*,*)  'Enter  title  of  file  containing  parameters' 

READ  (*,100)  im 

OPOKl,FIIf=TnL) 

READd.lOO)  TnT£ 

READ(1,110)  NPTS 
READ(1,110)  NRT 
READ(1,110)  Kd 
READ (1,110)  M 
READ  (1,110)  MLTAY 
READ (1,110)  NSTETPT 
READ(1,110)  NCAUS 
READ(1,100)  TTTU) 

READ(1,110)  lOffTS 
READ(1,110)  Kn 
READ(1,110)  INSTETPT 
CUDSE(l) 

GO  TO  215 
QiSE 
GO  TO  35 
BIDIF 

WRITE  (*,*)  '  ’ 


WRITE  (*,*)  'ISiter  title  of  file  containing  excitation  waveform' 


READ  (*,100)  TniD 
0PIN(8,FIl£=TniD) 

READ(8,100)  HOTES 

READ(8,U0)  N 

IF  (N  .CT.  1024)  IHEK 

VRITC  (*,*)  '»Afflber  of  points  in  data  file  exceer’s  the  dimensicn' 

WRITE  (»,*)  'of  the  array  used  in  the  program  to  store  the  file' 

STOP 

EMDIF 

CIiOSE{8) 

IF  ((N  .GE.  NHTS)  .AND.  DSCT)  THEN 

NDf>TS=N 

00  TO  215 

ENDIF 

NDPTS=N 

9  WRITE  (*,*)  'Enter  estimated  order  of  waveform' 

IF  (DSET)  THEN 

MAXMJHflPTSHM 

IF  (MAXIMUM  .or.  M-Kd-1)  MAXD«HH(d-l 

IF  (MAXIMUM  .or.  NDPTS-INSTRTTTHCn^l)  THEN 

MAXIMUM=NI»TrS-INSTTmT'-Kn-M+l 

ENDIF 

ELSE 

MAXIMUM=66 

ENDIF 

IF  (MAXIMUM  .EC.  1)  THEN 

WRITE  (*,*)  "The  estimated  order  of  the  waveform  can  only  be  1' 
IF  (DSET)  GO  TO  215 
GO  TO  10 
ELSE 

IF  (DSET)  THEN 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,’ 

ENDIF 

45  WRITE  (*,*)  'the  order  nay  range  from  1' 

WRITE  (*,*)  '  to’.MAXDtJM 

READ  (*,*)  Kn 

IF  (Kn  .GE.  1  .AND.  Kn  .LE.  MAXIMUM)  THEN 
IF  (DSET)  00  TO  215 
00  TO  10 
ENDIF 

WRITE  (*,*)  'Ehter  estimated  order  ajrain' 

WRITE  (*,*)  '  ' 

GO  TO  45 
ENDIF 

IF  (DSET)  00  TO  215 

10  INSTTm>T=l 

55  IF  (INSTHTPTtKn+M-1  .OT.  MOTS)  THEN 
INSTHTOT=INSTHTOT-1 
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ELSE 

IlCmPT=INSTOTm-l 

GO  TO  55 

EMDIF 

MSTOT=INSTOm' 


It  iJSffIkTFi  .iQ.  1)  THTO 

VRirc  (*,*)  "Hie  first  point  can  only  be  1' 

00  TO  215 
ELSE 

WRITE  {*,*)  'Diter  first  point  in  wavefora  file  to  be  processed* 
65  WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WRITE  {*,*)  'the  starting  point  may  range  frcm  1' 

WRITE  (*,*)  '  to',HSTRT 

READ  (*,*)  DCTRTPT 

IF  (INSnmT  .GE.  1  .AND.  INSTETPT  .LE.  MSTRT)  THfM 
IT  (DSET)  GO  TO  215 
GO  TO  1 
HJDIF 

WRITE  (*,*)  'filter  starting  point  again' 

WRITE  (*,*)  '  ' 

GO  TO  65 
n©IF 

IF  (DSET)  GO  TO  215 

1  IF  (.NOT.  DSET)  NUFnJB=.TEUE. 

IF  (.NOT.  DSET)  NOTRTPTV=1 

WRITE  (*,*)  'Biter  title  of  data  file  to  be  read' 

READ  (*,100)  TTITiE 
0P!N(12,FII£:Tni£) 

READ(12,100)  HEAKR 
READ (12, 110)  NPTS 
IF  (NPTS  .GT.  1024)  THEN 

WRITE  (*,*)  'Nimber  of  points  in  data  file  exceeds  the  dimension' 

WRITE  (*,*)  'of  the  array  used  in  the  program  to  store  the  file' 

STOP 

ENDIF 

ajOSE(12) 

IF  (NUFnz;  THEN 
GO  TO  3 

ELSEIF  (NSTRTm(Kd+«-l)*r*LTAY  .li).  NPTS)  THEN 

00  TO  215 

HjSE 

00  TO  6 

ENDIF 


3  IF  (NUFILE)  THEN 
MAXMJ»=69-Kn-1 
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IF  (MAXIMUM  .CT.  NPTS-69)  MAXDOHJPTS-69 
MI»=2 

IF  (MIN  MAXIMUM)  TOIN 
K(H(IN 

KRTTC  (*,*)  'Given  the  other  parameters  chosen  thus  far, ' 
WRITC  (*,*)  'Kd  must  be  ' .ruti 
GO  TO  4 
BiDIF 

VRITC  (*,*)  'Biter  Kd,  >=  the  estimated  order  of  the  system  ' 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

75  WRITE  (*,*)  'Kd  may  range  from', KIN 

WRITE  (*,*)  '  to', MAXIMUM 

READ  (*,*)  Kd 

IF  (Kd  .GE.  MIN  .AND.  Kd  .LE.  MAXIMUM)  GO  TO  4 
GO  TO  75 

BLSEIF  (DSET) 

MAXIMU«t-Kn-l 

IF  (MAXIMUM  .GT.  NPTS^l)  MAXIMUM=NPTS41 

KI»=2 

N=MAXIMUM 

85  IF  (NSTTm>T+(N+M-l)*raiiTAY  .LE.  NPTS)  THEN 
MAXDUHI 

IF  (MIN  .R3.  MAXIMUM)  TTON 

Kd=MIN 

GO  TO  215 

ELSEIF  (MAXIMUM  .LT.  MIN)  THEN 
liLTAY=l 

IF  (l+(2+M-l)*EeiTAY  .l£.  NPTS)  THEN 
Kd-2 

GO  TO  135 
ENDIF 

WRITE  {*,*)  'Error.  Kd  must  be  less  than  2' 

Kd=2 

GO  TO  215 
ENDIF 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

95  WRITE  (*,*)  'Kd  may  range  from  ',MIN 

WRITE  {*,*)  '  to' .MAXIMUM 

WRITE  {*,*)  'Eliter  Kd' 

READ  {*,*)  Kd 

IF  (Kd  .GE.  MIN  .AND.  Kd  .!£.  MAXIMUM)  GO  TO  215 
GO  TO  95 
ELSE 
EHJ-1 
00  TO  85 
ENDIF 
ENDIF 
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r 


I 

t 


C  Determine  M 

4  IF  (NUFILE)  TOHJ 

WRTIi:  {*,*)  'Biter  M,  the  row  dimension  of  the  data  natrix' 
IF  (.NOT.  DSEn*  .»0).  LONS)  THEN 
WRTTC  (*,*)  '  ' 

NRITC  (*,*)  'Note:  Kd+M  points  in  '.title 
HRITC  (*,*)  '  will  be  processed  ' 

WRITC  (*,*)  '  ' 

ENDIF 

105  HRTTE  (*,*)  'M  may  range  froo'.Kd 
IF  (NPTS-Kd  .CT.  69)  THQJ 

WRITC  (*,*)  '  to  69' 

WRITE  (*,*)  '  to', NPTS-Kd 

HffilF 

READ  (*,*)  M 
IF  (M  .CT.  69)  *raQ( 

WRITE  (*,*)  'M  must  also  be  less  than  70' 

GO  It)  105 

ELSEIF  (M  .LT.  Kd)  TEEN 

WRITE  (*,*)  'M  must  be  greater  than  or  equal  to  Kd,  Kd=  ',Kd 
GO  TO  105 

ELSEIF  {Kd4M  .GT.  NPTS)  THEN 

WRITE  (*,*)  'Rdm  must  be  less  than  or  equal  to'  ,NPTS, ' , ' 
WRITE  {‘,*)  'the  nunber  of  data  points  in'.TnUB 
WRITE  (*,*)  '  ' 

GO  ID  105 
ENDIF 

C  Begin  part  for  data  already  set 
ELSE 
ENCd 

115  IF  {NSTETPTf(Kd+«-l)*raLTAY  .l£.  NPTS)  THEN 
N=Ntl 
GO  TO  115 
PICT! 
fH»-l 
ENDIF 

IF  (N  .EC.  Kd)  TEEN 

WRITE  (*,*)  'M  must  equal', Kd 

N=Kd 

GO  TO  215 

ENDIF 

MAXDUHi 

IF  (MAXDfM  .CT.  69)  NAXIHUH=«9 
IF  (Kd+Kn+1  .IQ.  MAXDUl)  TEEN 
N=Kd+ftH-l 
00  TO  215 

ELSEIF  (KdWCiH-l  .OT.  MAXMM)  TEEN 
WRITE  (*,*)  'Kd  must  be  reduced' 
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GO  TO  3 
ELSE 

IlIM=KiWto+l 

EMDIF 

IE  (MIN  .LT.  ItoHCd+1)  MIN=ito+Kd+l 
125  VRmi  (*,*)  'M  may  range  f ran' , KIN 
TOTTE  (*,‘)  ’  to'.MMOJO! 

NRITC  (*,*)  'Ehter  M’ 

READ  (*,*)  M 

IF  (M  .GE.  MIN  .AND.  M  .IZ,  MAXIMUM)  GO  TO  215 

GO  TO  125 

QIDIF 


135 

5 

145 


155 


Determine  EELTAY 

IF  {.NOT.  NUFIIZ)  00  TO  215 

N=1 

IF  (NSTOTPT+N*(Kd4M-l)  .IZ.  NPTS)  TWM 
N=*H1 
00  TO  145 
ELSE 
N^l 
INDIF 

IF  (N  .EJ3,  1)  THEM 

WRITC  (*,*)  'Given  the  other  parameters  chosen  thus  far,’ 
WRITE  (*,*)  'Spacing  can  only  be  1' 

EELTAY=1 


IF  (NUFUZ)  THEN 
GO  TO  165 
ELSE 

GO  TO  215 
ENDIF 
ENDIF 

IF  (.NOT.  DSET  .AND.  LCNG)  THEN 

WRITE  (*,*)  'Ehter  spacing  between  the  ',Kd+M 

WRITE  (*,*)  'data  points  of  '.TTnZ 

WRITE  (*,*)  'to  be  processed  ' 

WRITE  (*,*)  ’  ' 

WRITE  (*,*) 

WRITE  {*.*) 

WRITE  (*,*) 

WRITE  (*,*)  ’  ’ 

OZE 

WRITE  {*,*)  'Biter  spacing  ' 

WRITE  (*,*)  '  ' 

ENDIF 

WRITE  (*,*)  'Spacing  may  range  fron 
WRITE  (*,*)  '  to',N 

READ  (*,*)  ULTAY 

IF  (EHiTAY  .a.  1  .AND.  KLTAY  .IZ.  N)  TEEN 
IF  (MJFHZ)  TEEN 


'If,  for  example,  one  is  chosen,  then  ’,Kd+M 
'consecutive  points  in  ',TIT1Z 
'will  be  processed  ' 
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GO  TO  165 
ELSE 

00  TO  215 
EMDIF 
OiSE 

00  TO  155 
EMDIF 

165  WRrni  (*,*)  'Do  you  wish  to  adjust  eigenvalues?  (y/n) ' 

RERD  (*,150)  YN 

IF  (YN  .EC.  'N'  .CR.  YN  .E®.  'n') 

IF  (NUFn£)  00  TO  6 

00  TO  215 

ENDIF 

IF  (YN  .NE.  'Y'  .Mffi.  YN  .NE.  'y')  GO  TO  165 
2  WRITE  (*,*)  'Discard  or  ccn^iensate  eigenvalues?  (d/c) ' 

RERD  (*,150)  DC 

IF  (DC  .EC.  'D'  .OR.  DC  .EC-  'd')  TTSl 

NRT=Kd 

00  TO  175 

EKDIF 

IF  (DC  .NE.  'C  .RND.  DC  .NE.  'c')  00  TO  2 

WRITE  (*,*)  'ELter  estijnate  of  the  actual  order  of  the  system' 

WRITE  (*,*)  '  ' 

IF  (LCNG)  TTIEN 

WRITE  (*,*)  "Ihis  estimate  will  be  used  to  determine  the  ' 
WRITE  (*,*)  'number  of  eigeivalues  cco^aensated  or  discarded  ’ 
EHDIF 

175  WRITE  (*,*)  'the  estimate  may  range  from  2' 

WRITE  (*,*)  '  to',Kd+Kn+l 

REM  (*,*)  NRT 

IF  (NRT  .CT.  KiMn+l  .OR.  NRT  .LT.  2)  THEH 
GO  TO  175 

ELSEIF  (.NOT.  NUFILE)  THEH 

GO  TO  215 

ENDIF 

6  1CTRTPT=1 

185  IF  (reTRTm(Kd+«-l)*I*LTRY  .!£.  NPTS)  THEN 
NSTRm^=NS7ETPT+l 
00  TO  185 
ELSE 

NSTRTTT=NSTRTPT-1 

EKDIF 

IF  (NSnmT  .EC.  1)  TEEJJ 

WRITE  (*,*)  'Given  the  other  parameters  chosen  thus  far,' 

WRITE  (*,*)  'the  starting  point  for  processing  the  data' 

WRITE  (*,*)  'must  be  the  first  point  in  the  data  file' 

00  TO  215 
EMDIF 
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WRITE  (*,*)  'Biter  desired  starting  point  in  data  file' 

IF  (.NOT.  DSBT  .AND.  LONG)  HON 

WRITE  {*,*)  '1  indicates  the  first  point  in  the  data  file  ' 
BIDIF 

wmE  (*,*)  '  ' 

WRITE  {*,*)  'Giver  the  other  parameters  chosen  thus  far,' 
195  WRITE  (*,*)  'the  starting  point  may  range  frco  1' 

WRITE  (*,*)  '  to',N?nmT 

READ  (*,*)  N 

IF  (N  .GE.  1  .AND.  N  .!£.  NSTRTTT)  TMW 

Nsram^ 

QiSE 

WRITE  (*,*)  'Biter  starting  point  again' 

WRITE  {*,*)  '  ' 

00  TO  195 
B0)IF 

IF  (.NOT.  NUnii:)  GO  TO  215 

7  IF  (DSET)  T«N 

IF  (NCAUS  .B2.  1)  THEN 

NCAUS=2 

00  TO  215 

ELSE 

NCAUS=1 

00  TO  215 

BffilF 

ENDIF 

WRITE  (*,*)  'Do  you  want  the  data  matrix  arrangement  to  be' 
WRITE  (*,’‘)  '  ’ 

WRITE  (*,*)  '1.  Causal' 

WRITE  (*,*)  '2.  Non-causal' 

WRITE  (*,*)  '  ' 

205  WRITE  (*,*)  'Please  enter  1  or  2  ' 

READ  (*,*)  NCAUS 
IF  (NCAUS  .EC.  1)  THEN 
CAUSAL=.TRUE. 

ELSEIF  (NCAUS  .EQ.  2)  THEN 
CAUSA]>.FALSE. 

ELSE 

00  TO  205 
INDIF 
00  TO  215 


12  WRITE  (*,*)  'Biter  title  of  file  to  contain  parameters' 
REM  (*,100)  im 
CPEN(l,FIli>=TnL) 

WRITE(1,100)  Tmi: 

WRITE (1,110)  NPTS 
WRITE(1,110)  NRT 
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I 


MRrn:{i,uo)  Kd 

WRITE  (1,110)  M 
WRITE  (1,110)  EQTRY 
WRITE  (1,110)  fCTRTTT 
WRITE  (1,110)  NCWJS 
WRITE(1,100)  TTITi) 
WRITE(1,110)  NDPTS 
WRITE  (1,110)  Kn 
WRITE  (1,110)  INSTWET 
aOSE(l) 

IF  (DSET)  GO  It)  215 

15  IF  (DSET)  TB!N 
CL0SE(2) 

CliOSEO) 

CALL  SUBPLT(NCVFRLAY) 
INDIF 


215  DSEr=.TRUE. 

MJFIU^  •  FALSE « 

WRnE(*,'‘)  '  ■ 

WRITE(*,*)  '1.  Data  file  to  be  processed  ',T 

+mE 

WRITE(*,*)  '  Number  of  data  points  in  data  file  ’,NPTS 

WRnE(*,*)  '2.  Estimated  order  of  the  system  ',NRT 

WRnE(*,*)  '3.  Kd,  the  number  of  columns  in  the  data  matrix’, Kd 
WRITE(*,‘')  '4.  M,  the  nimiber  of  rows  in  the  data  matrix’, M 
WRITE(*,*)  ’5.  Spacing  between  data  points  being  processed  ’,I®L‘IA 
+Y 

WRnE(*,*)  ’6,  First  point  in  the  data  file  to  be  processed’ ,NSTRT 
+PT 

WRnE(*,*)  ’  Last  point  in  the  data  file  to  be  processed’ ,NSTRT 
+PT+Kdf«-1 

IF  (NCAUS  .B2.  1)  TWN 

WRnE(*,*)  ’7.  Data  matrix  arrangemsit  for  processing  CA 

■KJSAL  ’ 

ELSE 

WRITE (*,*)  '7.  Data  matrix  arrangement  for  processing  Ndf-CA 
■KJSAL  ’ 

E»DIF 

WRnE(*,*)  ’ 


WRnE(*,*)  ’8,  File  containing  excitation  waveform  ’,T 

■HTTi) 

WRriE(*,*)  ’  Number  of  data  points  in  above  file  ’, NDPTS 

WRnE(*,*)  ’9.  Estimated  order  of  the  waveform  ’,Kn 

WRITE(*,*)  ’10.  First  point  in  the  file  to  be  ’ 

WRnE(*,»)  ’  input  into  the  data  matrix  ’,INSTR 

■HTT 
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wRrn:(*,*) 

t  t 

MRrn:(»,*) 

'll. 

«Rrn:(*,*) 

'12, 

iiRrn:{*,*) 

'13. 

+  file' 

WRrra:(*,*) 

'14. 

WRrn:(*,*) 

'15. 

'16, 

■Kig' 

Begin  processing  using  above  settings' 

Store  parameters  1-10  in  a  file' 

Retrieve  parameters  1-10  frcm  a  previously  created 

Reset  overlays' 

Re-plot  overlays' 

Bid  this  session  of  Cadzcw-Solcmcn  signal  process! 


TOrn:(*,*)  '  ' 

WOTEC*,*)  'filter  an  integer  from  1  to  16  to  make  changes  as  often 


+  as  you  desire' 

READ  (»,*)  NMQJU 

IF  (KMENU  .LT.  1  .OR.  WENU  .GT.  16)  THIH 
WtrnK*,*)  'filter  an  integer  from  1  to  16' 
00  TO  225 
ENDIF 


GO  10  (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12,13, 14, 15, 16), WffiMJ 

0PEN(12,FIIi>=Trn£) 

READ (12, 100)  REAKR 
READ(12,110)  NPTS 
READ (12, 120)  XQ 
READ(12,120)  XQ 
DO  235  1=1,  NPTS 
READ(12,120)  D(I) 

OCMTINtJE 

aiDSE(12) 

0PIN(8,nii>=TraD) 

READ(6,100)  HEADfS 
READ(8,iio)  Nirrs 
READ(8,120)  XQ 
READ(8,120)  XQ 
DO  245  1=1, NOTTS 
READ(8,120)  DK(I) 
oarriNUE 
CLDSE(8) 

KdPLT=Kd 

TOrrE(*,*)  'enter  title  of  file  to  contain  real  part  of  poles' 

RnD(*,100)  TmiF 

0PI«(2,nil!=TnUR) 


WRrrE(*,*) 'enter  title  of  file  to  contain  imaginary  part  of  poles 

READ(*,100)  Trn£I 

OPOJ(3,nii;^TTLEI) 


WRrn:(10,130)  (KdPLT) 

wirre(io,ioo)  TruiR 

WRITE (10, 100)  TTIlZr 
130  P0RMM'(I2) 

IIH!AX(M,Kd+Kn+l) 

100  PCRIftT(A) 
no  FORMAT  (15) 

120  roRMAT(E12.6) 

150  FaaiAT(A) 


DO  255  I=l,Kd4« 

(I)  =D  ( (I-l)  *I»LTAY4NSTRTfT) 

255  OCNTMJE 

265  DO  285  1=1,  M 

DO  275  J=l,Kd-H(n+l 
A(I,J)=Py(I+J) 

IF  (J  .<3:.  Kd+1)  A(I,J)=Dx(I-K3+INSnmT-2-Kd) 

275  OCWTINUE 
285  CCOTINUE 

B(l)=Dy(l) 

DO  295  1=2, M 
B(I)=A(I-1,1) 

295  OCrriNUE 

N=i(d+Kn+1 

C  Begin  singular  value  decxnposition 

CALL  SVD(MACHEP,M,N,W,A,W,MA'nj,U,MATV,V,rEKR,RVl) 

C  Errors  in  SVD? 

IF  (lERR  .GT.  0.0)  11®! 

WRITE  (*,*)  'Error  in  singular  value  nunber  ’,IERR,STCP 
EUDIF 

IF  (YN  .5).  'N')  C»  TO  385 

DO  305  I=l,Kd+Kn+l 
XP(I)=C.O 
305  OCKTIMUE 

C  Discard  or  ocnpensate  eigenvalues 

c  Order  singular  values 

XP(1)=W(1) 
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DO  335  I=2,Kd+Kn+l 
DO  325  J=1,I 

if  (W(I)  .GT.  XP{J))  THEN 
DO  315  K=I+1,J,-1 
315  XP(K)=XP(K-1) 

XP(j)=4Ki) 

GO  TO  335 
FMDIF 

325  CCOTINUE 

XP{I+1)^(I) 

335  aWITNUE 

C  XP(  )  new  contains  ordered  singular  values:  XP(1)  is  the  largest 

C  Discard  eigenvalues 

IF  (DC  .EQ.  'D')  THEN 
DO  345  J=Mm-l,Kd+{Cn+l 
345  W(J)=(0.0) 

ELSE 

C  Corpensste  eigenvalues 

AVG=0.0 

DO  355  J=4®Tfl,Kd-H(rH-l 
AVG=AVG+XP(J) *»2 
355  CCMTINUE 

IF  (Kd+Kn+l  .GT.  NRT)  AVG=:AVG/DBIJE (FLOAT (Kd+Kn+l-NRT) ) 

DO  375  J=l,Kd+Kn+l 
DO  365  K=l,Kd-HCn+l 
IF  (  V(J)  .B3.  XP(K)  )  THEN 
IF  (  K  .CT.  NRT  )  ITOl 
W(J)=0.0 
ELSE 

W(J)=DS<3n'(DABS(  W(J)*W(J)-AVG)) 

ENDIF 
GO  TO  375 
ENDIF 
365  CXWTINUE 
375  OinTNUE 
ENDIF 


385  DO  405  1=1, « 

DO  395  J=1,M 
irr(I,J)  =  (U(J,I)) 

395  OCWriNUE 
405  OCWTINUE 

C  Form  SI01A+  (Kd+Kn+1  x  M) 
DO  425  I=l,Kd+Kn+l 
DO  415  J=1,M 
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SI(M(I,J)=0.0 

IF  (I  .B2.  J  -AND.  W(J)  .HE.  0.0)  TWU 
SI(J1A(I,J)=1.0dO/W{J) 

ELSE 

SIGMA(I,J)=O.QDO 

QJDIT 

415  OCWriNUE 
425  OCNTINUE 

C  Fboa  SIGMA  (M  x  Kd+Kn+1) 

DO  445  1=1,  M 
DO  435  J=l, Kd+Kn+1 
SIG{I,J)=0.0 

IF  (I  .B3.  J)  SIG(I,J)=W{J) 

435  OCWTNUE 
445  CCNTUWE 

C  V=Kd+K^^^lxKld+KI^^l,SI<M+<=Kd+Kr^^lxM,VS=4M 

CAUi  MXMIIL(V, SIGMA, Ka+Kn+l,Kd+Kn+l,M, VS) 

C  VS=*i+Kn+b«,^rM^xM,AIN\NCd+Kr^^lx^ 

CALL  MXMLJL(VS,irr,Kd+Kn+l,M,M,AINV) 

C  Calculate  matrix  multiplication  of  AINV  x  B,  where 

C  AIN\N<d+Ktrt-lxM,B=Mxl,XP=Kd+Kn+lxl 

CAa  MXMUL(AINV,B,Kfl+Knr+l,M,L,XP) 

C  Conpute  autoregressive  coefficients  from  prediction  coefficients 

IF  (XP(Kd)  .B3.  0.0)  THfJ) 

WRITE  (*,*)  'E3?R0R,  avoiding  division  by  zero' 

STIF 

ELSE 

B(Kd)=1.0d0/XP(Kd) 

ENDIF 

DO  455  1=2, Kd 
B(I-l)=-B(Kd)*XP(Kd-i+l) 

455  OQOTINUE 

DO  465  i=l,Kd 
X(I)=-B(Kd-I+l) 

IF  (NCAUS  .FQ.  1)  X(I)=-XP(Kd-I+l) 

465  OCWnNUE 

X(Kd+l)=1.0 

C  Conpute  the  nxits  of  the  polynomial  in  z 

CAIA  PCI[i?T(X,0Cf’,KD,R0Crre,R0Cm,IIR) 

IF  (HR  .NE.  0)  WRITE  (*,+)  'fUROR  with  PCU?r,  nR=',m,ST0P 


138 


DO  475  1=1, Kd 
WRITE(2,120)  ROCTR(I) 

VRrrE(3,120)  ROOTI(I) 

S(I)=DCMPLX(ROCTO(I)  ,ROCm(I) ) 

475  aamnuE 

MAGfCLO 
DO  485  1=1, Kd 

IT  (a»BS(S{I))  .GE.  l.ODO)  MAGKti=4ftGra/+l 
485  CCKTINUE 

WRnE(*,*)  ’#  of  poles  Kith  magnitude  <=  1' ,Kd-MAGPOIi 
VRITE  (*,*)  'HCT  ANY  KEY  TO  CCmNUE* 

READ  (*,100)  HEADS? 

C  Plot  poles 

NCVD?1AY=N0VB?LAY+1 

ai0SE(2) 

ciiOSEii; 

CALL  SlBPLT(NOVB?LAY) 

J=0 

K=0 

DO  495  1=1, Kd 

IF  (CDABS^S(I')  .LT.  1.0)  ’11©) 

WRITE  (*,*)  S(I) ,CDABS(S{I)) 

J=J+1 

K=K+1 

BJDIF 

IF  (J  .B?.  20)  THDJ 

WRITE  (*,*)  'HIT  ANY  KEY  TO  COOTINUE' 

READ  (*,100)  HEADS? 

J=0 

BJDIF 

495  OCWITNUE 

WRITE  (*,*)  'Poles  with  magnitude  less  than  cne  ',K 
WRITE  (*,*)  'HIT  ANY  KEY  TO  CCKTINUE' 

READ  (*,100)  HEAIS? 

00  TO  215 

16  STOP 
BJD 

'Z 
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APPENDIX  C.  MATRIX  MULTIPLICATION 


SBROUITNE  M3MJL(A,B,RA,CA,CB,AB) 

BmXSR  RA,CA,CB 

RZMi*8  An0,70),B(70,70),AB{70,70) 

C  Calculates  matrix  multiplication  of  A  x  B=AB,  where 

C  A=«AxCA,B=CAxCB,AB==RAxCB 

DO  30  1=1, RA 
DO  20  >1=1,  CB 
AB(I,J)=O.0 
DO  10  K=1,CA 

AB(I,J)=AB(I,J)+A(I,K)*B(K,J) 

10  OOriTNUE 

20  ccwmwE 

30  CCWITNUE 

RETURN 

Dffi 
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APPENDIX  D.  GRAPHICS  ROUTINE 


SIBROUITNE  aBPLT(NWnUAy) 


c 

C  HS-FCRTRAN  Program  using  "Grafmatic"  Lilarary  Subroutines. 

C  Plots  a  Solid  Line  and  Optional  Overlay  Plot  for  Cmgarison. 

C  Written  by  M.A.  Morgan  with  Latest  Update  August  1989. 

C 

C  Default  Printer  is  "IBM  Graphics"  (e.g.  Edison,  (Sddata,  IBM) 

C  With  Plot  Rotated  90  degrees  Fran  the  Vertical.  "GrafPlus.Ccm" 

C  May  be  Run  to  Rotate  Plot  Upright  cn  Paper  and  to  Use  a  Variety 

C  of  Intact  Printers.  "GrafLaser.Com"  May  be  Run  to  Use  a  Laser 
C  Printer.  See  Graf Plus/Laser  Manual  FTom  Jewell  Technology. 

C 

c 

CHARACTER*!  YN,¥Nl,DUM,YN2,SYMBCXi,BEiL,FEE2),FFYN 
CHARACnR*4  LINE 
CHARACTER*?  SVMB 

CHARACTTR*16  LTTT,CTrr,FTlAME,TriUR,TITLEI 
CHARACnR*64  Tmi:,HCCFY 
REAL  CRTR{70),CRn(70),NRTR(70),NRTI{70) 

INTBa5R*2  N,JRCW,JC^AJSYKl,ISYM2,ITYPEl,ITyPE2,NSCRN 
IWrEX3R*2  CyAN,GREEN,WHm,YELU»,RIB,BLACK,BLUE,NTW3 
INm3R*2  JROWl,JROW2,JCm,JCaj2,CROSS,KdPLT,I 
INIEX2R*2  PURPLE,  RUST 
EXTERNAL  XFUN,YFUNP,YFTMJ 
LINE=' —  ’ 
wHrrE=7 
GRIH»=10 
CYAN=11 
YELLCIW=14 
RH>=12 
BLACK=0 
BLUE=1 
NTHD=2 
PWPIi>=5 
RUST=€ 

■m=CHAR{7) 

FEH)=CHAR(12) 

C  Gear  Screen  and  Put  Up  Introduction  -  on  Blue  BacAgound  for  EJGA 
C  Only;  Another  Background  Color  is  Possible  by  Changing  "BdilT' 

C  in  the  Calls  to  gPRBG  and  QCWSCN. 

CALL  QSMOra:(NTVO) 

CALL  Qf>RBG(0,BIiJE) 
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CALL  QCVSCN(BUJE) 

WRITE!*,*)  BEL 
NS=1 

NSCRN  =16 
lTyPE2=<) 

Calling  GKAmATIC  Routines  and  Plotting  FI  Solid  Line  Oraph 
lTyPEl=0 
1SYM1=-1 
NDOTS1=0 
JR0W1=1 
JR0W2=350 
J0CL1=  75 
J0CL2=  565 
XmN^l.2 
XMAX=1.2 
YMIN=-1.2 
VMAX=1.20 
YOVm=1.115 
XDRG=0.0 
VORG=0.0 
XST=-1.1 
XFIN=1.1 
YST=-1.1 
Yn»=l.l 

CAIL  QSMXE  (NSCRN) 

CALL  (yiOT(JOCM,JCXMi2,JRCWl,JROW2,XMIN,XMAX,YMIN,Ym,»DRG,ra^ 
+l,TOVnOC,l,5) 

CAIL  QSEITJP(ND0TS1, CYAN, ISYm, RED) 

IF(XFIN-XST  .LE.  9.0)  XMAJ0R=O.6 
IF(XFIN-XST  .l£.  6.0)  XMAJOR=0.4 
IF(XFIN-XST  .LE.  3.3)  XMAJOR=0.2 
rF{XFIN-XST  .GE.  9.0)  X10tK)R=(XnN-XST) /lO.O 

MI1WN3 

LABEL=1 

NEBC=2 

CAIL  QmiS(XST,XFIN,XMIWrcR,KDai,LA^ 
y>aX)R=XMAJCR 

CAIL  CranS(YST,YriN,YMAJOR,KINaR,LABEL,NITC 
Plot  unit  circle 
A=^■1.0 
B=1.0 

CAIL  QCURV(XFl«,YFl«P,A,B) 

CALL  QCURV(XFW,YFU«,A,B) 

IF  (NCWEWAY-l  .LT.  1)  THIN 
IF  {NDVI»LAY-1  .RJ.  0)  TOUJ 
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WRITE  (*,3)  NOVnOAY-l 
ELSE 

WRITC  (*,3)  NZH» 

ENDIF 

HiSEIF  (NCNfHOAY-l  .GT.  1)  TOTO 

WRITE  (*,3)  NOVIRIAY-1 

ELSE 

WRTIT:  (*,4)  NDVIRIAY-1 
ENDIF 

3  FORMAT  (13,'  OVERLAYS') 

4  ET*MAT  (13,'  OVERLAY  ') 

REWIND  (10) 


DO  20  I=l,NOVERLAY 
READ  (10,110)  KdPLT 
READ  (10,100)  TnUR 
READ  (10,100)  TITlzr 
QPEN(2,ni£=TrniR) 

OPEN  ( 3 ,  Fni>=Trn£r ) 
l«d=KdPLT 
DO  27  J=l, KdPLT 
READ  (2,120)  NRTR(J) 

READ  (3,120)  NRTI(J) 

IF  (DSQRT(M?rR(J)**2+NRn(J)**2)  .GT.  1.1)  TTON 

NKd=t«d-l 

NRTO(J)=0.0 

NRn(J)=0.0 

ENDIF 

27  ocrmNUE 
PURPLE=5 
RUST=€ 

WHITE=7 

GREEN=10 

CYAI^ll 

YEIiliOW=14 

RED=12 

BLUE>=1 

IF  (I  .B3.  1)  THEN 
(ALL  QSETUP(NDOrSl,CYAN,ISYMl,RED) 

ELSEIF  (I  .EC-  2)  TOEN 
CAIi  QSETOP(NDOI'Sl,(nfAN,ISYMl, GREEN) 

ELSCTF  (I  .EC.  3)  THEN 
CALL  QSETUP(NDOrSl,CYAN,ISYMl,YEUiCIW) 

ELSEITF  (I  .EC-  4)  THEN 
(m  QSEIUP(NDOTS1,CYAN,ISYIO.,BLUE) 

ELSEIF  (I  .EC-  5)  THEN 
CALL  QSErnjP(NDOTSl,CYAN,ISYMl,WHrTE) 

ELSEIF  (I  .EC-  6)  THEN 
CAa  QSETUF;NDOrSl,CYAN,ISYm,FURPI£) 
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giSEIF  (I  .DC.  7)  raUf 
cm  QSErup(Nix7rsi,cv»N,isnfl,RusT) 

ELSE 

cm  QSEIT)P(NDCfrSl,CY»N,ISYIIl,Rn)) 

EKDIF 

cm  c?im(mPEi,KdPLT,Nm,NRn) 

20  CXWriMJE 

READ  (*,100)  DOM 
GO  TO  40 

HC3CPy=' HARDCOPY — >  HfTSK  P  OR  p' 

CALL  Cpr!n’(30,HC0PY,RID,25,l) 
cm  QCM0V(55,1) 

Hccpy=' 

cm  QpT3n’{40,H0CPY, BLACK, 25,1) 

IF  (DUM  .HE.  'P'  .AND.  DUM  .NE.  ’p’)  GO  TO  40 

cm  ^«CRN 

OPJN  {1,FILE='PRN') 

WRITE  (1,160)  FOB 
100  roRMAT(A) 

110  P3?MAT(I2) 

120  roRMAT(E12.6) 

160  TORMATC  ’,A,\) 

40  OOmNUE 

cm  QSMa*:(Nrwo) 
cm  (yRBG(0,BUJE) 
cm  CPVSCN(BLUE) 

WRTIEC*,*)  NKd, 'points  were  plotted’ 

RETURN 

END 

REAL  FUNCnCfJ  XFUN(T) 

XFUN=T 

RETURN 

END 

REAL  FTJKCnCN  YEWJPd) 

YFUNP=SQRT(1.0-T*T) 

RETURN 

END 

REAL  FlWCnCN  YF(f«(T) 

YFUI*=-SQRT(1.0-T*T) 

RFIURN 

END 
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